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 !is_null(Ops_[numRowBlocks_ * j + i]);
235}
236
237
238template<class Scalar>
240 const int i, const int j
241 ) const
242{
243 if (!blockExists(i, j)) return false;
244 return Ops_[numRowBlocks_*j+i].isConst();
245}
246
247
248template<class Scalar>
251{
252 if (!blockExists(i, j)) return Teuchos::null;
253 return Ops_[numRowBlocks_*j+i].getNonconstObj();
254}
255
256
257template<class Scalar>
259DefaultBlockedLinearOp<Scalar>::getBlock(const int i, const int j) const
260{
261 if (!blockExists(i, j)) return Teuchos::null;
262 return Ops_[numRowBlocks_*j+i];
263}
264
265
266// Overridden from LinearOpBase
267
268
269template<class Scalar>
272{
273 return productRange_;
274}
275
276
277template<class Scalar>
280{
281 return productDomain_;
282}
283
284
285template<class Scalar>
288{
289 return Teuchos::null; // ToDo: Implement this when needed!
290}
291
292
293// Overridden from Teuchos::Describable
294
295
296template<class Scalar>
298{
299 assertBlockFillIsActive(false);
300 std::ostringstream oss;
301 oss
303 << "numRowBlocks="<<numRowBlocks_
304 << ",numColBlocks="<<numColBlocks_
305 << "}";
306 return oss.str();
307}
308
309
310template<class Scalar>
312 Teuchos::FancyOStream &out_arg
313 ,const Teuchos::EVerbosityLevel verbLevel
314 ) const
315{
316 using Teuchos::rcpFromRef;
318 using Teuchos::OSTab;
319 assertBlockFillIsActive(false);
320 RCP<FancyOStream> out = rcpFromRef(out_arg);
321 OSTab tab1(out);
322 switch(verbLevel) {
325 *out << this->description() << std::endl;
326 break;
330 {
331 *out
333 << "rangeDim=" << this->range()->dim()
334 << ",domainDim=" << this->domain()->dim()
335 << ",numRowBlocks=" << numRowBlocks_
336 << ",numColBlocks=" << numColBlocks_
337 << "}\n";
338 OSTab tab2(out);
339 *out
340 << "Constituent LinearOpBase objects for M = [ Op[0,0] ..."
341 << " ; ... ; ... Op[numRowBlocks-1,numColBlocks-1] ]:\n";
342 tab2.incrTab();
343 for( int i = 0; i < numRowBlocks_; ++i ) {
344 for( int j = 0; j < numColBlocks_; ++j ) {
345 *out << "Op["<<i<<","<<j<<"] = ";
347 block_i_j = getBlock(i,j);
348 if(block_i_j.get())
349 *out << Teuchos::describe(*getBlock(i,j),verbLevel);
350 else
351 *out << "NULL\n";
352 }
353 }
354 break;
355 }
356 default:
357 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
358 }
359}
360
361
362// protected
363
364
365// Overridden from LinearOpBase
366
367
368template<class Scalar>
370{
371 bool supported = true;
372 for( int i = 0; i < numRowBlocks_; ++i ) {
373 for( int j = 0; j < numColBlocks_; ++j ) {
375 block_i_j = getBlock(i,j);
376 if( block_i_j.get() && !Thyra::opSupported(*block_i_j,M_trans) )
377 supported = false;
378 }
379 }
380 return supported;
381}
382
383
384template<class Scalar>
386 const EOpTransp M_trans,
387 const MultiVectorBase<Scalar> &X_in,
388 const Ptr<MultiVectorBase<Scalar> > &Y_inout,
389 const Scalar alpha,
390 const Scalar beta
391 ) const
392{
393
394 using Teuchos::rcpFromRef;
396 typedef RCP<MultiVectorBase<Scalar> > MultiVectorPtr;
397 typedef RCP<const MultiVectorBase<Scalar> > ConstMultiVectorPtr;
398 typedef RCP<const LinearOpBase<Scalar> > ConstLinearOpPtr;
399
400#ifdef TEUCHOS_DEBUG
402 "DefaultBlockedLinearOp<Scalar>::apply(...)", *this, M_trans, X_in, &*Y_inout
403 );
404#endif // TEUCHOS_DEBUG
405
406 const bool
407 struct_transp = (real_trans(M_trans)!=NOTRANS); // Structural transpose?
408 const int
409 opNumRowBlocks = ( !struct_transp ? numRowBlocks_ : numColBlocks_ ),
410 opNumColBlocks = ( !struct_transp ? numColBlocks_ : numRowBlocks_ );
411
412 //
413 // Y = alpha * op(M) * X + beta*Y
414 //
415 // =>
416 //
417 // Y[i] = beta+Y[i] + sum(alpha*op(Op)[i,j]*X[j],j=0...opNumColBlocks-1)
418 //
419 // , for i=0...opNumRowBlocks-1
420 //
421
423 defaultProductRange_op = ( real_trans(M_trans)==NOTRANS
424 ? defaultProductRange_ : defaultProductDomain_ ),
425 defaultProductDomain_op = ( real_trans(M_trans)==NOTRANS
426 ? defaultProductDomain_ : defaultProductRange_ );
427
429 X = castOrCreateSingleBlockProductMultiVector<Scalar>(
430 defaultProductDomain_op, rcpFromRef(X_in));
432 Y = nonconstCastOrCreateSingleBlockProductMultiVector<Scalar>(
433 defaultProductRange_op, rcpFromPtr(Y_inout));
434
435 for( int i = 0; i < opNumRowBlocks; ++i ) {
436 MultiVectorPtr Y_i = Y->getNonconstMultiVectorBlock(i);
437 for( int j = 0; j < opNumColBlocks; ++j ) {
438 ConstLinearOpPtr
439 Op_i_j = ( !struct_transp ? getBlock(i,j) : getBlock(j,i) );
440 ConstMultiVectorPtr
441 X_j = X->getMultiVectorBlock(j);
442 if(j==0) {
443 if (nonnull(Op_i_j))
444 Thyra::apply(*Op_i_j, M_trans,* X_j, Y_i.ptr(), alpha, beta);
445 else
446 scale(beta, Y_i.ptr());
447 }
448 else {
449 if (nonnull(Op_i_j))
450 Thyra::apply(*Op_i_j, M_trans, *X_j, Y_i.ptr(), alpha, ST::one());
451 }
452 }
453 }
454}
455
456// Overridden from RowStatLinearOpBase
457
458template<class Scalar>
459bool
461 const RowStatLinearOpBaseUtils::ERowStat row_stat) const
462{
463 using Teuchos::rcpFromRef;
464 using Teuchos::rcp_dynamic_cast;
465 using Teuchos::dyn_cast;
466 //typedef Teuchos::ScalarTraits<Scalar> ST; // unused
467 typedef RowStatLinearOpBase<Scalar> RowStatOp;
468 typedef RCP<const LinearOpBase<Scalar> > ConstLinearOpPtr;
469 typedef const LinearOpBase<Scalar> ConstLinearOp;
470
471 // Figure out what needs to be check for each sub block to compute
472 // the require row stat
473 RowStatLinearOpBaseUtils::ERowStat
474 subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
475 subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
476 switch (row_stat) {
477 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
478 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
479 subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM;
480 subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
481 break;
482 case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
483 case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
484 subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
485 subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM;
486 break;
487 default:
489 }
490
491 // check each sub block for compatibility (null means zero,
492 // automatically compatible)
493 for( int i = 0; i < numRowBlocks_; ++i ) {
494 for( int j = 0; j < numColBlocks_; ++j ) {
495 ConstLinearOpPtr Op_i_j = getBlock(i,j);
496
497 if (nonnull(Op_i_j)) {
498 // pull out adjoint, and scaling
499 ConstLinearOp * Orig_i_j = 0;
500 Scalar scalar = 1.0;
501 EOpTransp transp = NOTRANS;
502 Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
503
504 const RowStatOp & row_stat_op = Teuchos::dyn_cast<const RowStatOp>(*Orig_i_j);
505
506 // sub block must also support the required row stat operation
507 RowStatLinearOpBaseUtils::ERowStat stat = subblk_stat;
508 if(transp==NOTRANS || transp==CONJ)
509 stat = subblk_stat;
510 else if(transp==TRANS || transp==CONJTRANS)
511 stat = subblk_trans_stat;
512
513 if(!row_stat_op.rowStatIsSupported(stat))
514 return false;
515 }
516 // else: sub block is null, "0" is good enough
517 }
518 }
519
520 return true;
521}
522
523template<class Scalar>
524void
526 const RowStatLinearOpBaseUtils::ERowStat row_stat,
527 const Teuchos::Ptr<VectorBase< Scalar> > &rowStatVec) const
528{
529 using Teuchos::rcpFromRef;
530 using Teuchos::rcpFromPtr;
531 using Teuchos::rcp_dynamic_cast;
532 using Teuchos::dyn_cast;
534 typedef RowStatLinearOpBase<Scalar> RowStatOp;
535 typedef RCP<const LinearOpBase<Scalar> > ConstLinearOpPtr;
536 typedef const LinearOpBase<Scalar> ConstLinearOp;
537 typedef RCP<VectorBase<Scalar> > VectorPtr;
538
539 // Figure out what needs to be check for each sub block to compute
540 // the require row stat
541 RowStatLinearOpBaseUtils::ERowStat
542 subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
543 subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
544 switch (row_stat) {
545 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
546 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
547 subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM;
548 subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
549 break;
550 case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
551 case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
552 subblk_stat = RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM;
553 subblk_trans_stat = RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM;
554 break;
555 default:
557 }
558
560 Y = rcp_dynamic_cast<ProductVectorBase<Scalar> >(rcpFromPtr(rowStatVec));
561
562 // using sub block row stat capability interrogate each for
563 // its constribution
564 for( int i = 0; i < numRowBlocks_; ++i ) {
565 VectorPtr blk_vec = Y->getNonconstVectorBlock(i);
566 put_scalar (ST::zero (), blk_vec.ptr ()); // clear vector just in case
567
568 for( int j = 0; j < numColBlocks_; ++j ) {
569 ConstLinearOpPtr Op_i_j = getBlock(i,j);
570
571 VectorPtr tmp_vec = createMember(Op_i_j->range());
572
573 put_scalar (ST::zero (), tmp_vec.ptr ()); // clear vector just in case
574
575 if (nonnull(Op_i_j)) {
576 // pull out adjoint, and scaling
577 ConstLinearOp * Orig_i_j = 0;
578 Scalar scalar = 1.0;
579 EOpTransp transp = NOTRANS;
580 Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
581
582 const RowStatOp & row_stat_op = Teuchos::dyn_cast<const RowStatOp>(*Orig_i_j);
583
584 // sub block must also support the required row stat operation
585 if(transp==NOTRANS || transp==CONJ)
586 row_stat_op.getRowStat(subblk_stat,tmp_vec.ptr());
587 else if(transp==TRANS || transp==CONJTRANS)
588 row_stat_op.getRowStat(subblk_trans_stat,tmp_vec.ptr());
589 else
590 { TEUCHOS_ASSERT(false); }
591
592 // some the temporary into the block vector
593 Vp_V(blk_vec.ptr(),*tmp_vec);
594 }
595 }
596 }
597
598 // do post processing as needed (take the inverse currently
599 switch (row_stat) {
600 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
601 case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
602 reciprocal(*rowStatVec,rowStatVec.ptr());
603 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
604 case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
605 break;
606 default:
608 }
609}
610
611// Overridden from ScaledLinearOpBase
612
613template<class Scalar>
614bool
616{
617 using Teuchos::rcp_dynamic_cast;
618 using Teuchos::dyn_cast;
619 typedef RCP<const LinearOpBase<Scalar> > LinearOpPtr;
620 typedef const LinearOpBase<Scalar> LinearOp;
621 typedef const ScaledLinearOpBase<Scalar> ScaledOp;
622
623 bool supported = true;
624 for( int i = 0; i < numRowBlocks_; ++i ) {
625 for( int j = 0; j < numColBlocks_; ++j ) {
626 LinearOpPtr Op_i_j = getBlock(i,j);
627
628 EOpTransp transp = NOTRANS;
629
630 if (nonnull(Op_i_j)) {
631 // pull out adjoint, and scaling
632 LinearOp * Orig_i_j = 0;
633 Scalar scalar = 1.0;
634 Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
635
636 ScaledOp & scaled_op = Teuchos::dyn_cast<ScaledOp>(*Orig_i_j);
637
638 if(transp==NOTRANS || transp==CONJ)
639 supported &= scaled_op.supportsScaleLeft();
640 else if(transp==TRANS || transp==CONJTRANS)
641 supported &= scaled_op.supportsScaleRight();
642 }
643 }
644 }
645
646 return supported;
647}
648
649template<class Scalar>
650bool
652{
653 using Teuchos::rcp_dynamic_cast;
654 using Teuchos::dyn_cast;
655 typedef RCP<const LinearOpBase<Scalar> > LinearOpPtr;
656 typedef const LinearOpBase<Scalar> LinearOp;
657 typedef const ScaledLinearOpBase<Scalar> ScaledOp;
658
659 bool supported = true;
660 for( int i = 0; i < numRowBlocks_; ++i ) {
661 for( int j = 0; j < numColBlocks_; ++j ) {
662 LinearOpPtr Op_i_j = getBlock(i,j);
663
664 if (nonnull(Op_i_j)) {
665 // pull out adjoint, and scaling
666 LinearOp * Orig_i_j = 0;
667 Scalar scalar = 1.0;
668 EOpTransp transp = NOTRANS;
669 Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
670
671 ScaledOp & scaled_op = Teuchos::dyn_cast<ScaledOp>(*Orig_i_j);
672
673 if(transp==NOTRANS || transp==CONJ)
674 supported &= scaled_op.supportsScaleRight();
675 else if(transp==TRANS || transp==CONJTRANS)
676 supported &= scaled_op.supportsScaleLeft();
677 }
678 }
679 }
680
681 return supported;
682}
683
684template<class Scalar>
685void
687 const VectorBase< Scalar > &row_scaling)
688{
689 using Teuchos::dyn_cast;
690 using Teuchos::rcp_dynamic_cast;
691 typedef ScaledLinearOpBase<Scalar> ScaledOp;
692 typedef RCP<LinearOpBase<Scalar> > LinearOpPtr;
693 typedef RCP<const VectorBase<Scalar> > VectorPtr;
694
696 Y = dyn_cast<const ProductVectorBase<Scalar> >(row_scaling);
697
698 // using sub block row stat capability interrogate each for
699 // its constribution
700 for( int i = 0; i < numRowBlocks_; ++i ) {
701 VectorPtr blk_vec = Y.getVectorBlock(i);
702
703 for( int j = 0; j < numColBlocks_; ++j ) {
704 LinearOpPtr Op_i_j = getNonconstBlock(i,j);
705
706 if (nonnull(Op_i_j)) {
707 // pull out adjoint, and scaling
708 LinearOpPtr Orig_i_j;
709 EOpTransp transp = NOTRANS;
710
711 {
713 = rcp_dynamic_cast<ScaledAdjointLinearOpBase<Scalar> >(Op_i_j);
714 if(nonnull(saOp)) {
715 transp = saOp->overallTransp();
716 Orig_i_j = saOp->getNonconstOrigOp();
717 }
718 else
719 Orig_i_j = Op_i_j; // stick with the original
720 }
721
722 RCP<ScaledOp> scaled_op = rcp_dynamic_cast<ScaledOp>(Orig_i_j);
723
724 // sub block must support a row stat operation
725 TEUCHOS_ASSERT(scaled_op!=Teuchos::null); // guranteed by compatibility check
726
727 // sub block must also support the required row stat operation
728 if(transp==NOTRANS || transp==CONJ)
729 scaled_op->scaleLeft(*blk_vec);
730 else if(transp==TRANS || transp==CONJTRANS)
731 scaled_op->scaleRight(*blk_vec);
732 }
733 }
734 }
735}
736
737template<class Scalar>
738void
740 const VectorBase< Scalar > &col_scaling)
741{
742 using Teuchos::dyn_cast;
743 using Teuchos::rcp_dynamic_cast;
744 typedef ScaledLinearOpBase<Scalar> ScaledOp;
745 typedef RCP<LinearOpBase<Scalar> > LinearOpPtr;
746 typedef RCP<const VectorBase<Scalar> > VectorPtr;
747
749 Y = dyn_cast<const ProductVectorBase<Scalar> >(col_scaling);
750
751 // using sub block row stat capability interrogate each for
752 // its constribution
753 for( int j = 0; j < numColBlocks_; ++j ) {
754 VectorPtr blk_vec = Y.getVectorBlock(j);
755
756 for( int i = 0; i < numRowBlocks_; ++i ) {
757 LinearOpPtr Op_i_j = getNonconstBlock(i,j);
758
759 if (nonnull(Op_i_j)) {
760 // pull out adjoint, and scaling
761 LinearOpPtr Orig_i_j;
762 EOpTransp transp = NOTRANS;
763
764 {
766 = rcp_dynamic_cast<ScaledAdjointLinearOpBase<Scalar> >(Op_i_j);
767 if(nonnull(saOp)) {
768 transp = saOp->overallTransp();
769 Orig_i_j = saOp->getNonconstOrigOp();
770 }
771 else
772 Orig_i_j = Op_i_j; // stick with the original
773 }
774
775 RCP<ScaledOp> scaled_op = rcp_dynamic_cast<ScaledOp>(Orig_i_j);
776
777 // sub block must support a row stat operation
778 TEUCHOS_ASSERT(scaled_op!=Teuchos::null); // guranteed by compatibility check
779
780 // sub block must also support the required row stat operation
781 if(transp==NOTRANS || transp==CONJ)
782 scaled_op->scaleRight(*blk_vec);
783 else if(transp==TRANS || transp==CONJTRANS)
784 scaled_op->scaleLeft(*blk_vec);
785 }
786 }
787 }
788}
789
790// private
791
792
793template<class Scalar>
795 const int numRowBlocks, const int numColBlocks
796 )
797{
798 numRowBlocks_ = numRowBlocks;
799 numColBlocks_ = numColBlocks;
800 Ops_.resize(numRowBlocks_*numColBlocks_);
801 if (is_null(productRange_)) {
802 rangeBlocks_.resize(numRowBlocks);
803 domainBlocks_.resize(numColBlocks);
804 }
805 blockFillIsActive_ = true;
806}
807
808
809template<class Scalar>
810void DefaultBlockedLinearOp<Scalar>::assertBlockFillIsActive(
811 bool wantedValue
812 ) const
813{
814#ifdef TEUCHOS_DEBUG
815 TEUCHOS_TEST_FOR_EXCEPT(!(blockFillIsActive_==wantedValue));
816#else
817 (void)wantedValue;
818#endif
819}
820
821
822template<class Scalar>
823void DefaultBlockedLinearOp<Scalar>::assertBlockRowCol(
824 const int i, const int j
825 ) const
826{
827#ifdef TEUCHOS_DEBUG
829 !( 0 <= i ), std::logic_error
830 ,"Error, i="<<i<<" is invalid!"
831 );
833 !( 0 <= j ), std::logic_error
834 ,"Error, j="<<j<<" is invalid!"
835 );
836 // Only validate upper range if the number of row and column blocks is
837 // fixed!
838 if(Ops_.size()) {
840 !( 0 <= i && i < numRowBlocks_ ), std::logic_error
841 ,"Error, i="<<i<<" does not fall in the range [0,"<<numRowBlocks_-1<<"]!"
842 );
844 !( 0 <= j && j < numColBlocks_ ), std::logic_error
845 ,"Error, j="<<j<<" does not fall in the range [0,"<<numColBlocks_-1<<"]!"
846 );
847 }
848#else
849 (void)i;
850 (void)j;
851#endif
852}
853
854
855template<class Scalar>
856void DefaultBlockedLinearOp<Scalar>::setBlockSpaces(
857 const int i, const int j, const LinearOpBase<Scalar> &block
858 )
859{
860 using Teuchos::toString;
861 assertBlockFillIsActive(true);
862 assertBlockRowCol(i,j);
863
864 // Validate that if the vector space block is already set that it is
865 // compatible with the block that is being set.
866 if( i < numRowBlocks_ && j < numColBlocks_ ) {
867#ifdef TEUCHOS_DEBUG
868 RCP<const VectorSpaceBase<Scalar> >
869 rangeBlock = (
870 productRange_.get()
871 ? productRange_->getBlock(i)
872 : rangeBlocks_[i]
873 ),
874 domainBlock = (
875 productDomain_.get()
876 ? productDomain_->getBlock(j)
877 : domainBlocks_[j]
878 );
879 if(rangeBlock.get()) {
881 "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
882 "Adding block: " + block.description(),
883 *rangeBlock,("(*productRange->getBlock("+toString(i)+"))"),
884 *block.range(),("(*block["+toString(i)+","+toString(j)+"].range())")
885 );
886 }
887 if(domainBlock.get()) {
889 "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
890 "Adding block: " + block.description(),
891 *domainBlock,("(*productDomain->getBlock("+toString(j)+"))"),
892 *block.domain(),("(*block["+toString(i)+","+toString(j)+"].domain())")
893 );
894 }
895#endif // TEUCHOS_DEBUG
896 }
897
898 // Add spaces missing range and domain space blocks if we are doing a
899 // flexible fill (otherwise these loops will not be executed)
900 for( int k = numRowBlocks_; k <= i; ++k )
901 rangeBlocks_.push_back(Teuchos::null);
902 for( int k = numColBlocks_; k <= j; ++k )
903 domainBlocks_.push_back(Teuchos::null);
904
905 // Set the incoming range and domain blocks if not already set
906 if(!productRange_.get()) {
907 if(!rangeBlocks_[i].get())
908 rangeBlocks_[i] = block.range().assert_not_null();
909 if(!domainBlocks_[j].get()) {
910 domainBlocks_[j] = block.domain().assert_not_null();
911 }
912 }
913
914 // Update the current number of row and columns blocks if doing a flexible
915 // fill.
916 if(!Ops_.size()) {
917 numRowBlocks_ = rangeBlocks_.size();
918 numColBlocks_ = domainBlocks_.size();
919 }
920
921}
922
923
924template<class Scalar>
925template<class LinearOpType>
926void DefaultBlockedLinearOp<Scalar>::setBlockImpl(
927 const int i, const int j,
928 const RCP<LinearOpType> &block
929 )
930{
931 setBlockSpaces(i, j, *block);
932 if (Ops_.size()) {
933 // We are doing a fill with a fixed number of row and column blocks so we
934 // can just set this.
935 Ops_[numRowBlocks_*j+i] = block;
936 }
937 else {
938 // We are doing a flexible fill so add the block to the stack of blocks or
939 // replace a block that already exists.
940 bool foundBlock = false;
941 for( unsigned int k = 0; k < Ops_stack_.size(); ++k ) {
942 BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
943 if( block_i_j.i == i && block_i_j.j == j ) {
944 block_i_j.block = block;
945 foundBlock = true;
946 break;
947 }
948 }
949 if(!foundBlock)
950 Ops_stack_.push_back(BlockEntry<Scalar>(i,j,block));
951 }
952}
953
954
955template<class Scalar>
956void DefaultBlockedLinearOp<Scalar>::adjustBlockSpaces()
957{
958
959#ifdef TEUCHOS_DEBUG
960 TEUCHOS_ASSERT_INEQUALITY(Ops_.size(), !=, 0);
961#endif
962
963 //
964 // Loop through the rows and columns looking for rows with mixed
965 // single-space range and/or domain spaces on operators and set the single
966 // spaces as the block space if it exists.
967 //
968 // NOTE: Once we get here, we can safely assume that all of the operators
969 // are compatible w.r.t. their spaces so if there are rows and/or columns
970 // with mixed product and single vector spaces that we can just pick the
971 // single vector space for the whole row and/or column.
972 //
973
974 // Adjust blocks in the range space
975 for (int i = 0; i < numRowBlocks_; ++i) {
976 for (int j = 0; j < numColBlocks_; ++j) {
977 const RCP<const LinearOpBase<Scalar> >
978 op_i_j = Ops_[numRowBlocks_*j+i];
979 if (is_null(op_i_j))
980 continue;
981 const RCP<const VectorSpaceBase<Scalar> > range_i_j = op_i_j->range();
982 if (is_null(productVectorSpaceBase<Scalar>(range_i_j, false))) {
983 rangeBlocks_[i] = range_i_j;
984 break;
985 }
986 }
987 }
988
989 // Adjust blocks in the domain space
990 for (int j = 0; j < numColBlocks_; ++j) {
991 for (int i = 0; i < numRowBlocks_; ++i) {
992 const RCP<const LinearOpBase<Scalar> >
993 op_i_j = Ops_[numRowBlocks_*j+i];
994 if (is_null(op_i_j))
995 continue;
996 const RCP<const VectorSpaceBase<Scalar> >
997 domain_i_j = op_i_j->domain();
998 if (is_null(productVectorSpaceBase<Scalar>(domain_i_j, false))) {
999 domainBlocks_[j] = domain_i_j;
1000 break;
1001 }
1002 }
1003 }
1004
1005}
1006
1007
1008} // namespace Thyra
1009
1010
1011template<class Scalar>
1013Thyra::defaultBlockedLinearOp()
1014{
1015 return Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
1016}
1017
1018
1019template<class Scalar>
1021Thyra::block1x1(
1022 const RCP<const LinearOpBase<Scalar> > &A00,
1023 const std::string &label
1024 )
1025{
1026 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1027 M = defaultBlockedLinearOp<Scalar>();
1028 M->beginBlockFill(1,1);
1029 M->setBlock(0, 0, A00);
1030 M->endBlockFill();
1031 if (label.length())
1032 M->setObjectLabel(label);
1033 return M;
1034}
1035
1036
1037template<class Scalar>
1039Thyra::block1x2(
1040 const RCP<const LinearOpBase<Scalar> > &A00,
1041 const RCP<const LinearOpBase<Scalar> > &A01,
1042 const std::string &label
1043 )
1044{
1045 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1046 M = defaultBlockedLinearOp<Scalar>();
1047 M->beginBlockFill(1,2);
1048 if(A00.get()) M->setBlock(0,0,A00);
1049 if(A01.get()) M->setBlock(0,1,A01);
1050 M->endBlockFill();
1051 if (label.length())
1052 M->setObjectLabel(label);
1053 return M;
1054}
1055
1056
1057template<class Scalar>
1059Thyra::block2x1(
1060 const RCP<const LinearOpBase<Scalar> > &A00,
1061 const RCP<const LinearOpBase<Scalar> > &A10,
1062 const std::string &label
1063 )
1064{
1065 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1066 M = defaultBlockedLinearOp<Scalar>();
1067 M->beginBlockFill(2,1);
1068 if(A00.get()) M->setBlock(0,0,A00);
1069 if(A10.get()) M->setBlock(1,0,A10);
1070 M->endBlockFill();
1071 if (label.length())
1072 M->setObjectLabel(label);
1073 return M;
1074}
1075
1076
1077template<class Scalar>
1079Thyra::block2x2(
1080 const RCP<const LinearOpBase<Scalar> > &A00,
1081 const RCP<const LinearOpBase<Scalar> > &A01,
1082 const RCP<const LinearOpBase<Scalar> > &A10,
1083 const RCP<const LinearOpBase<Scalar> > &A11,
1084 const std::string &label
1085 )
1086{
1087 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1088 M = defaultBlockedLinearOp<Scalar>();
1089 M->beginBlockFill(2,2);
1090 if(A00.get()) M->setBlock(0,0,A00);
1091 if(A01.get()) M->setBlock(0,1,A01);
1092 if(A10.get()) M->setBlock(1,0,A10);
1093 if(A11.get()) M->setBlock(1,1,A11);
1094 M->endBlockFill();
1095 if (label.length())
1096 M->setObjectLabel(label);
1097 return M;
1098}
1099
1100
1101template<class Scalar>
1103Thyra::nonconstBlock1x1(
1104 const RCP<LinearOpBase<Scalar> > &A00,
1105 const std::string &label
1106 )
1107{
1108 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1109 M = defaultBlockedLinearOp<Scalar>();
1110 M->beginBlockFill(1, 1);
1111 M->setNonconstBlock(0, 0, A00);
1112 M->endBlockFill();
1113 if (label.length())
1114 M->setObjectLabel(label);
1115 return M;
1116}
1117
1118
1119template<class Scalar>
1121Thyra::nonconstBlock1x2(
1122 const RCP<LinearOpBase<Scalar> > &A00,
1123 const RCP<LinearOpBase<Scalar> > &A01,
1124 const std::string &label
1125 )
1126{
1127 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1128 M = defaultBlockedLinearOp<Scalar>();
1129 M->beginBlockFill(1,2);
1130 if(A00.get()) M->setNonconstBlock(0,0,A00);
1131 if(A01.get()) M->setNonconstBlock(0,1,A01);
1132 M->endBlockFill();
1133 if (label.length())
1134 M->setObjectLabel(label);
1135 return M;
1136}
1137
1138
1139template<class Scalar>
1141Thyra::nonconstBlock2x1(
1142 const RCP<LinearOpBase<Scalar> > &A00,
1143 const RCP<LinearOpBase<Scalar> > &A10,
1144 const std::string &label
1145 )
1146{
1147 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1148 M = defaultBlockedLinearOp<Scalar>();
1149 M->beginBlockFill(2,1);
1150 if(A00.get()) M->setNonconstBlock(0,0,A00);
1151 if(A10.get()) M->setNonconstBlock(1,0,A10);
1152 M->endBlockFill();
1153 if (label.length())
1154 M->setObjectLabel(label);
1155 return M;
1156}
1157
1158
1159template<class Scalar>
1161Thyra::nonconstBlock2x2(
1162 const RCP<LinearOpBase<Scalar> > &A00,
1163 const RCP<LinearOpBase<Scalar> > &A01,
1164 const RCP<LinearOpBase<Scalar> > &A10,
1165 const RCP<LinearOpBase<Scalar> > &A11,
1166 const std::string &label
1167 )
1168{
1169 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1170 M = defaultBlockedLinearOp<Scalar>();
1171 M->beginBlockFill(2,2);
1172 if(A00.get()) M->setNonconstBlock(0,0,A00);
1173 if(A01.get()) M->setNonconstBlock(0,1,A01);
1174 if(A10.get()) M->setNonconstBlock(1,0,A10);
1175 if(A11.get()) M->setNonconstBlock(1,1,A11);
1176 M->endBlockFill();
1177 if (label.length())
1178 M->setObjectLabel(label);
1179 return M;
1180}
1181
1182
1183//
1184// Explicit instantiation macro
1185//
1186// Must be expanded from within the Thyra namespace!
1187//
1188
1189
1190#define THYRA_DEFAULT_BLOCKED_LINEAR_OP_INSTANT(SCALAR) \
1191 \
1192 template class DefaultBlockedLinearOp<SCALAR >; \
1193 \
1194 template RCP<DefaultBlockedLinearOp<SCALAR > > \
1195 defaultBlockedLinearOp<SCALAR >(); \
1196 \
1197 template RCP<const LinearOpBase<SCALAR > > \
1198 block1x1( \
1199 const RCP<const LinearOpBase<SCALAR > > &A00, \
1200 const std::string &label \
1201 ); \
1202 \
1203 template RCP<const LinearOpBase<SCALAR > > \
1204 block1x2( \
1205 const RCP<const LinearOpBase<SCALAR > > &A00, \
1206 const RCP<const LinearOpBase<SCALAR > > &A01, \
1207 const std::string &label \
1208 ); \
1209 \
1210 template RCP<const LinearOpBase<SCALAR > > \
1211 block2x1( \
1212 const RCP<const LinearOpBase<SCALAR > > &A00, \
1213 const RCP<const LinearOpBase<SCALAR > > &A10, \
1214 const std::string &label \
1215 ); \
1216 \
1217 template RCP<const LinearOpBase<SCALAR > > \
1218 block2x2( \
1219 const RCP<const LinearOpBase<SCALAR > > &A00, \
1220 const RCP<const LinearOpBase<SCALAR > > &A01, \
1221 const RCP<const LinearOpBase<SCALAR > > &A10, \
1222 const RCP<const LinearOpBase<SCALAR > > &A11, \
1223 const std::string &label \
1224 ); \
1225 \
1226 template RCP<LinearOpBase<SCALAR > > \
1227 nonconstBlock1x1( \
1228 const RCP<LinearOpBase<SCALAR > > &A00, \
1229 const std::string &label \
1230 ); \
1231 \
1232 template RCP<LinearOpBase<SCALAR > > \
1233 nonconstBlock1x2( \
1234 const RCP<LinearOpBase<SCALAR > > &A00, \
1235 const RCP<LinearOpBase<SCALAR > > &A01, \
1236 const std::string &label \
1237 ); \
1238 \
1239 template RCP<LinearOpBase<SCALAR > > \
1240 nonconstBlock2x1( \
1241 const RCP<LinearOpBase<SCALAR > > &A00, \
1242 const RCP<LinearOpBase<SCALAR > > &A10, \
1243 const std::string &label \
1244 ); \
1245 \
1246 template RCP<LinearOpBase<SCALAR > > \
1247 nonconstBlock2x2( \
1248 const RCP<LinearOpBase<SCALAR > > &A00, \
1249 const RCP<LinearOpBase<SCALAR > > &A01, \
1250 const RCP<LinearOpBase<SCALAR > > &A10, \
1251 const RCP<LinearOpBase<SCALAR > > &A11, \
1252 const std::string &label \
1253 );
1254
1255
1256#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)