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>
31DefaultBlockedLinearOp<Scalar>::DefaultBlockedLinearOp()
32 :numRowBlocks_(0), numColBlocks_(0), blockFillIsActive_(false)
33{}
34
35
36// Overridden from PhysicallyBlockedLinearOpBase
37
38
39template<class Scalar>
40void DefaultBlockedLinearOp<Scalar>::beginBlockFill()
41{
42 assertBlockFillIsActive(false);
43 uninitialize();
44 resetStorage(0,0);
45}
46
47
48template<class Scalar>
49void DefaultBlockedLinearOp<Scalar>::beginBlockFill(
50 const int numRowBlocks, const int numColBlocks
51 )
52{
53 assertBlockFillIsActive(false);
54 uninitialize();
55 resetStorage(numRowBlocks,numColBlocks);
56}
57
58
59template<class Scalar>
60void DefaultBlockedLinearOp<Scalar>::beginBlockFill(
61 const RCP<const ProductVectorSpaceBase<Scalar> > &new_productRange
62 ,const RCP<const ProductVectorSpaceBase<Scalar> > &new_productDomain
63 )
64{
65 using Teuchos::rcp_dynamic_cast;
66 assertBlockFillIsActive(false);
67 uninitialize();
68 productRange_ = new_productRange.assert_not_null();
69 productDomain_ = new_productDomain.assert_not_null();
70 defaultProductRange_ =
71 rcp_dynamic_cast<const DefaultProductVectorSpace<Scalar> >(productRange_);
72 defaultProductDomain_ =
73 rcp_dynamic_cast<const DefaultProductVectorSpace<Scalar> >(productDomain_);
74 // Note: the above spaces must be set before calling the next function!
75 resetStorage(productRange_->numBlocks(), productDomain_->numBlocks());
76}
77
78
79template<class Scalar>
80bool DefaultBlockedLinearOp<Scalar>::blockFillIsActive() const
81{
82 return blockFillIsActive_;
83}
84
85
86template<class Scalar>
87bool DefaultBlockedLinearOp<Scalar>::acceptsBlock(
88 const int i, const int j
89 ) const
90{
91 assertBlockFillIsActive(true);
92 assertBlockRowCol(i,j);
93 return true;
94}
95
96
97template<class Scalar>
98void DefaultBlockedLinearOp<Scalar>::setNonconstBlock(
99 const int i, const int j
100 ,const RCP<LinearOpBase<Scalar> > &block
101 )
102{
103 setBlockImpl(i, j, block);
104}
105
106
107template<class Scalar>
108void DefaultBlockedLinearOp<Scalar>::setBlock(
109 const int i, const int j
110 ,const RCP<const LinearOpBase<Scalar> > &block
111 )
112{
113 setBlockImpl(i, j, block);
114}
115
116
117template<class Scalar>
118void DefaultBlockedLinearOp<Scalar>::endBlockFill()
119{
120
121 using Teuchos::as;
122
123 assertBlockFillIsActive(true);
124
125 // 2009/05/06: rabartl: ToDo: When doing a flexible block fill
126 // (Ops_stack_.size() > 0), we need to assert that all of the block rows and
127 // columns have been filled in. I don't think we do that here.
128
129 // Get the number of block rows and columns
130 if (nonnull(productRange_)) {
131 numRowBlocks_ = productRange_->numBlocks();
132 numColBlocks_ = productDomain_->numBlocks();
133 }
134 else {
135 numRowBlocks_ = rangeBlocks_.size();
136 numColBlocks_ = domainBlocks_.size();
137 // NOTE: Above, whether doing a flexible fill or not, all of the blocks
138 // must be set in order to have a valid filled operator so this
139 // calculation should be correct.
140 }
141
142 // Assert that all of the block rows and columns have at least one entry if
143 // the spaces were not given up front.
144#ifdef TEUCHOS_DEBUG
145 if (is_null(productRange_)) {
146 for (int i = 0; i < numRowBlocks_; ++i) {
147 TEUCHOS_TEST_FOR_EXCEPTION(
148 !rangeBlocks_[i].get(), std::logic_error
149 ,"DefaultBlockedLinearOp<Scalar>::endBlockFill():"
150 " Error, no linear operator block for the i="<<i<<" block row was added"
151 " and we can not complete the block fill!"
152 );
153 }
154 for(int j = 0; j < numColBlocks_; ++j) {
155 TEUCHOS_TEST_FOR_EXCEPTION(
156 !domainBlocks_[j].get(), std::logic_error
157 ,"DefaultBlockedLinearOp<Scalar>::endBlockFill():"
158 " Error, no linear operator block for the j="
159 <<j<<" block column was added"
160 " and we can not complete the block fill!"
161 );
162 }
163 }
164#endif
165
166 // Insert the block LOB objects if doing a flexible fill.
167 if (Ops_stack_.size()) {
168 Ops_.resize(numRowBlocks_*numColBlocks_);
169 for ( int k = 0; k < as<int>(Ops_stack_.size()); ++k ) {
170 const BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
171 Ops_[numRowBlocks_*block_i_j.j + block_i_j.i] = block_i_j.block;
172 }
173 Ops_stack_.resize(0);
174 }
175
176 // Set the product range and domain spaces if not already set
177 if (is_null(productRange_)) {
178 adjustBlockSpaces();
179 defaultProductRange_ = productVectorSpace<Scalar>(rangeBlocks_());
180 defaultProductDomain_ = productVectorSpace<Scalar>(domainBlocks_());
181 productRange_ = defaultProductRange_;
182 productDomain_ = defaultProductDomain_;
183 }
184
185 rangeBlocks_.resize(0);
186 domainBlocks_.resize(0);
187
188 blockFillIsActive_ = false;
189
190}
191
192
193template<class Scalar>
194void DefaultBlockedLinearOp<Scalar>::uninitialize()
195{
196 productRange_ = Teuchos::null;
197 productDomain_ = Teuchos::null;
198 numRowBlocks_ = 0;
199 numColBlocks_ = 0;
200 Ops_.resize(0);
201 Ops_stack_.resize(0);
202 rangeBlocks_.resize(0);
203 domainBlocks_.resize(0);
204 blockFillIsActive_ = false;
205}
206
207
208// Overridden from BlockedLinearOpBase
209
210
211template<class Scalar>
212RCP<const ProductVectorSpaceBase<Scalar> >
213DefaultBlockedLinearOp<Scalar>::productRange() const
214{
215 return productRange_;
216}
217
218
219template<class Scalar>
220RCP<const ProductVectorSpaceBase<Scalar> >
221DefaultBlockedLinearOp<Scalar>::productDomain() const
222{
223 return productDomain_;
224}
225
226
227template<class Scalar>
228bool DefaultBlockedLinearOp<Scalar>::blockExists(
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>
239bool DefaultBlockedLinearOp<Scalar>::blockIsConst(
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>
249RCP<LinearOpBase<Scalar> >
250DefaultBlockedLinearOp<Scalar>::getNonconstBlock(const int i, const int j)
251{
252 if (!blockExists(i, j)) return Teuchos::null;
253 return Ops_[numRowBlocks_*j+i].getNonconstObj();
254}
255
256
257template<class Scalar>
258RCP<const LinearOpBase<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>
270RCP< const VectorSpaceBase<Scalar> >
271DefaultBlockedLinearOp<Scalar>::range() const
272{
273 return productRange_;
274}
275
276
277template<class Scalar>
278RCP< const VectorSpaceBase<Scalar> >
279DefaultBlockedLinearOp<Scalar>::domain() const
280{
281 return productDomain_;
282}
283
284
285template<class Scalar>
286RCP<const LinearOpBase<Scalar> >
287DefaultBlockedLinearOp<Scalar>::clone() const
288{
289 return Teuchos::null; // ToDo: Implement this when needed!
290}
291
292
293// Overridden from Teuchos::Describable
294
295
296template<class Scalar>
297std::string DefaultBlockedLinearOp<Scalar>::description() const
298{
299 assertBlockFillIsActive(false);
300 std::ostringstream oss;
301 oss
302 << Teuchos::Describable::description() << "{"
303 << "numRowBlocks="<<numRowBlocks_
304 << ",numColBlocks="<<numColBlocks_
305 << "}";
306 return oss.str();
307}
308
309
310template<class Scalar>
311void DefaultBlockedLinearOp<Scalar>::describe(
312 Teuchos::FancyOStream &out_arg
313 ,const Teuchos::EVerbosityLevel verbLevel
314 ) const
315{
316 using Teuchos::rcpFromRef;
317 using Teuchos::FancyOStream;
318 using Teuchos::OSTab;
319 assertBlockFillIsActive(false);
320 RCP<FancyOStream> out = rcpFromRef(out_arg);
321 OSTab tab1(out);
322 switch(verbLevel) {
323 case Teuchos::VERB_DEFAULT:
324 case Teuchos::VERB_LOW:
325 *out << this->description() << std::endl;
326 break;
327 case Teuchos::VERB_MEDIUM:
328 case Teuchos::VERB_HIGH:
329 case Teuchos::VERB_EXTREME:
330 {
331 *out
332 << Teuchos::Describable::description() << "{"
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<<"] = ";
346 RCP<const LinearOpBase<Scalar> >
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>
369bool DefaultBlockedLinearOp<Scalar>::opSupportedImpl(EOpTransp M_trans) const
370{
371 bool supported = true;
372 for( int i = 0; i < numRowBlocks_; ++i ) {
373 for( int j = 0; j < numColBlocks_; ++j ) {
374 RCP<const LinearOpBase<Scalar> >
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>
385void DefaultBlockedLinearOp<Scalar>::applyImpl(
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;
395 typedef Teuchos::ScalarTraits<Scalar> ST;
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
401 THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
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
422 const RCP<const DefaultProductVectorSpace<Scalar> >
423 defaultProductRange_op = ( real_trans(M_trans)==NOTRANS
424 ? defaultProductRange_ : defaultProductDomain_ ),
425 defaultProductDomain_op = ( real_trans(M_trans)==NOTRANS
426 ? defaultProductDomain_ : defaultProductRange_ );
427
428 const RCP<const ProductMultiVectorBase<Scalar> >
429 X = castOrCreateSingleBlockProductMultiVector<Scalar>(
430 defaultProductDomain_op, rcpFromRef(X_in));
431 const RCP<ProductMultiVectorBase<Scalar> >
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
460DefaultBlockedLinearOp<Scalar>::rowStatIsSupportedImpl(
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:
488 TEUCHOS_TEST_FOR_EXCEPT(true);
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
525DefaultBlockedLinearOp<Scalar>::getRowStatImpl(
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;
533 typedef Teuchos::ScalarTraits<Scalar> ST;
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:
556 TEUCHOS_TEST_FOR_EXCEPT(true);
557 }
558
559 const RCP<ProductVectorBase<Scalar> >
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 // a null block is zero and does not contribute
572 if (nonnull(Op_i_j)) {
573 VectorPtr tmp_vec = createMember(Op_i_j->range());
574
575 put_scalar (ST::zero (), tmp_vec.ptr ()); // clear vector just in case
576
577 // pull out adjoint, and scaling
578 ConstLinearOp * Orig_i_j = 0;
579 Scalar scalar = 1.0;
580 EOpTransp transp = NOTRANS;
581 Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
582
583 const RowStatOp & row_stat_op = Teuchos::dyn_cast<const RowStatOp>(*Orig_i_j);
584
585 // sub block must also support the required row stat operation
586 if(transp==NOTRANS || transp==CONJ)
587 row_stat_op.getRowStat(subblk_stat,tmp_vec.ptr());
588 else if(transp==TRANS || transp==CONJTRANS)
589 row_stat_op.getRowStat(subblk_trans_stat,tmp_vec.ptr());
590 else
591 { TEUCHOS_ASSERT(false); }
592
593 // some the temporary into the block vector
594 Vp_V(blk_vec.ptr(),*tmp_vec);
595 }
596 }
597 }
598
599 // do post processing as needed (take the inverse currently
600 switch (row_stat) {
601 case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
602 case RowStatLinearOpBaseUtils::ROW_STAT_INV_COL_SUM:
603 reciprocal(*rowStatVec,rowStatVec.ptr());
604 case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
605 case RowStatLinearOpBaseUtils::ROW_STAT_COL_SUM:
606 break;
607 default:
608 TEUCHOS_TEST_FOR_EXCEPT(true);
609 }
610}
611
612// Overridden from ScaledLinearOpBase
613
614template<class Scalar>
615bool
616DefaultBlockedLinearOp<Scalar>::supportsScaleLeftImpl() const
617{
618 using Teuchos::rcp_dynamic_cast;
619 using Teuchos::dyn_cast;
620 typedef RCP<const LinearOpBase<Scalar> > LinearOpPtr;
621 typedef const LinearOpBase<Scalar> LinearOp;
622 typedef const ScaledLinearOpBase<Scalar> ScaledOp;
623
624 bool supported = true;
625 for( int i = 0; i < numRowBlocks_; ++i ) {
626 for( int j = 0; j < numColBlocks_; ++j ) {
627 LinearOpPtr Op_i_j = getBlock(i,j);
628
629 EOpTransp transp = NOTRANS;
630
631 if (nonnull(Op_i_j)) {
632 // pull out adjoint, and scaling
633 LinearOp * Orig_i_j = 0;
634 Scalar scalar = 1.0;
635 Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
636
637 ScaledOp & scaled_op = Teuchos::dyn_cast<ScaledOp>(*Orig_i_j);
638
639 if(transp==NOTRANS || transp==CONJ)
640 supported &= scaled_op.supportsScaleLeft();
641 else if(transp==TRANS || transp==CONJTRANS)
642 supported &= scaled_op.supportsScaleRight();
643 }
644 }
645 }
646
647 return supported;
648}
649
650template<class Scalar>
651bool
652DefaultBlockedLinearOp<Scalar>::supportsScaleRightImpl() const
653{
654 using Teuchos::rcp_dynamic_cast;
655 using Teuchos::dyn_cast;
656 typedef RCP<const LinearOpBase<Scalar> > LinearOpPtr;
657 typedef const LinearOpBase<Scalar> LinearOp;
658 typedef const ScaledLinearOpBase<Scalar> ScaledOp;
659
660 bool supported = true;
661 for( int i = 0; i < numRowBlocks_; ++i ) {
662 for( int j = 0; j < numColBlocks_; ++j ) {
663 LinearOpPtr Op_i_j = getBlock(i,j);
664
665 if (nonnull(Op_i_j)) {
666 // pull out adjoint, and scaling
667 LinearOp * Orig_i_j = 0;
668 Scalar scalar = 1.0;
669 EOpTransp transp = NOTRANS;
670 Thyra::unwrap(*Op_i_j,&scalar,&transp,&Orig_i_j);
671
672 ScaledOp & scaled_op = Teuchos::dyn_cast<ScaledOp>(*Orig_i_j);
673
674 if(transp==NOTRANS || transp==CONJ)
675 supported &= scaled_op.supportsScaleRight();
676 else if(transp==TRANS || transp==CONJTRANS)
677 supported &= scaled_op.supportsScaleLeft();
678 }
679 }
680 }
681
682 return supported;
683}
684
685template<class Scalar>
686void
687DefaultBlockedLinearOp<Scalar>::scaleLeftImpl(
688 const VectorBase< Scalar > &row_scaling)
689{
690 using Teuchos::dyn_cast;
691 using Teuchos::rcp_dynamic_cast;
692 typedef ScaledLinearOpBase<Scalar> ScaledOp;
693 typedef RCP<LinearOpBase<Scalar> > LinearOpPtr;
694 typedef RCP<const VectorBase<Scalar> > VectorPtr;
695
696 const ProductVectorBase<Scalar> &
697 Y = dyn_cast<const ProductVectorBase<Scalar> >(row_scaling);
698
699 // using sub block row stat capability interrogate each for
700 // its constribution
701 for( int i = 0; i < numRowBlocks_; ++i ) {
702 VectorPtr blk_vec = Y.getVectorBlock(i);
703
704 for( int j = 0; j < numColBlocks_; ++j ) {
705 LinearOpPtr Op_i_j = getNonconstBlock(i,j);
706
707 if (nonnull(Op_i_j)) {
708 // pull out adjoint, and scaling
709 LinearOpPtr Orig_i_j;
710 EOpTransp transp = NOTRANS;
711
712 {
713 RCP<ScaledAdjointLinearOpBase<Scalar> > saOp
714 = rcp_dynamic_cast<ScaledAdjointLinearOpBase<Scalar> >(Op_i_j);
715 if(nonnull(saOp)) {
716 transp = saOp->overallTransp();
717 Orig_i_j = saOp->getNonconstOrigOp();
718 }
719 else
720 Orig_i_j = Op_i_j; // stick with the original
721 }
722
723 RCP<ScaledOp> scaled_op = rcp_dynamic_cast<ScaledOp>(Orig_i_j);
724
725 // sub block must support a row stat operation
726 TEUCHOS_ASSERT(scaled_op!=Teuchos::null); // guranteed by compatibility check
727
728 // sub block must also support the required row stat operation
729 if(transp==NOTRANS || transp==CONJ)
730 scaled_op->scaleLeft(*blk_vec);
731 else if(transp==TRANS || transp==CONJTRANS)
732 scaled_op->scaleRight(*blk_vec);
733 }
734 }
735 }
736}
737
738template<class Scalar>
739void
740DefaultBlockedLinearOp<Scalar>::scaleRightImpl(
741 const VectorBase< Scalar > &col_scaling)
742{
743 using Teuchos::dyn_cast;
744 using Teuchos::rcp_dynamic_cast;
745 typedef ScaledLinearOpBase<Scalar> ScaledOp;
746 typedef RCP<LinearOpBase<Scalar> > LinearOpPtr;
747 typedef RCP<const VectorBase<Scalar> > VectorPtr;
748
749 const ProductVectorBase<Scalar> &
750 Y = dyn_cast<const ProductVectorBase<Scalar> >(col_scaling);
751
752 // using sub block row stat capability interrogate each for
753 // its constribution
754 for( int j = 0; j < numColBlocks_; ++j ) {
755 VectorPtr blk_vec = Y.getVectorBlock(j);
756
757 for( int i = 0; i < numRowBlocks_; ++i ) {
758 LinearOpPtr Op_i_j = getNonconstBlock(i,j);
759
760 if (nonnull(Op_i_j)) {
761 // pull out adjoint, and scaling
762 LinearOpPtr Orig_i_j;
763 EOpTransp transp = NOTRANS;
764
765 {
766 RCP<ScaledAdjointLinearOpBase<Scalar> > saOp
767 = rcp_dynamic_cast<ScaledAdjointLinearOpBase<Scalar> >(Op_i_j);
768 if(nonnull(saOp)) {
769 transp = saOp->overallTransp();
770 Orig_i_j = saOp->getNonconstOrigOp();
771 }
772 else
773 Orig_i_j = Op_i_j; // stick with the original
774 }
775
776 RCP<ScaledOp> scaled_op = rcp_dynamic_cast<ScaledOp>(Orig_i_j);
777
778 // sub block must support a row stat operation
779 TEUCHOS_ASSERT(scaled_op!=Teuchos::null); // guranteed by compatibility check
780
781 // sub block must also support the required row stat operation
782 if(transp==NOTRANS || transp==CONJ)
783 scaled_op->scaleRight(*blk_vec);
784 else if(transp==TRANS || transp==CONJTRANS)
785 scaled_op->scaleLeft(*blk_vec);
786 }
787 }
788 }
789}
790
791// private
792
793
794template<class Scalar>
795void DefaultBlockedLinearOp<Scalar>::resetStorage(
796 const int numRowBlocks, const int numColBlocks
797 )
798{
799 numRowBlocks_ = numRowBlocks;
800 numColBlocks_ = numColBlocks;
801 Ops_.resize(numRowBlocks_*numColBlocks_);
802 if (is_null(productRange_)) {
803 rangeBlocks_.resize(numRowBlocks);
804 domainBlocks_.resize(numColBlocks);
805 }
806 blockFillIsActive_ = true;
807}
808
809
810template<class Scalar>
811void DefaultBlockedLinearOp<Scalar>::assertBlockFillIsActive(
812 bool wantedValue
813 ) const
814{
815#ifdef TEUCHOS_DEBUG
816 TEUCHOS_TEST_FOR_EXCEPT(!(blockFillIsActive_==wantedValue));
817#else
818 (void)wantedValue;
819#endif
820}
821
822
823template<class Scalar>
824void DefaultBlockedLinearOp<Scalar>::assertBlockRowCol(
825 const int i, const int j
826 ) const
827{
828#ifdef TEUCHOS_DEBUG
829 TEUCHOS_TEST_FOR_EXCEPTION(
830 !( 0 <= i ), std::logic_error
831 ,"Error, i="<<i<<" is invalid!"
832 );
833 TEUCHOS_TEST_FOR_EXCEPTION(
834 !( 0 <= j ), std::logic_error
835 ,"Error, j="<<j<<" is invalid!"
836 );
837 // Only validate upper range if the number of row and column blocks is
838 // fixed!
839 if(Ops_.size()) {
840 TEUCHOS_TEST_FOR_EXCEPTION(
841 !( 0 <= i && i < numRowBlocks_ ), std::logic_error
842 ,"Error, i="<<i<<" does not fall in the range [0,"<<numRowBlocks_-1<<"]!"
843 );
844 TEUCHOS_TEST_FOR_EXCEPTION(
845 !( 0 <= j && j < numColBlocks_ ), std::logic_error
846 ,"Error, j="<<j<<" does not fall in the range [0,"<<numColBlocks_-1<<"]!"
847 );
848 }
849#else
850 (void)i;
851 (void)j;
852#endif
853}
854
855
856template<class Scalar>
857void DefaultBlockedLinearOp<Scalar>::setBlockSpaces(
858 const int i, const int j, const LinearOpBase<Scalar> &block
859 )
860{
861 using Teuchos::toString;
862 assertBlockFillIsActive(true);
863 assertBlockRowCol(i,j);
864
865 // Validate that if the vector space block is already set that it is
866 // compatible with the block that is being set.
867 if( i < numRowBlocks_ && j < numColBlocks_ ) {
868#ifdef TEUCHOS_DEBUG
869 RCP<const VectorSpaceBase<Scalar> >
870 rangeBlock = (
871 productRange_.get()
872 ? productRange_->getBlock(i)
873 : rangeBlocks_[i]
874 ),
875 domainBlock = (
876 productDomain_.get()
877 ? productDomain_->getBlock(j)
878 : domainBlocks_[j]
879 );
880 if(rangeBlock.get()) {
881 THYRA_ASSERT_VEC_SPACES_NAMES(
882 "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
883 "Adding block: " + block.description(),
884 *rangeBlock,("(*productRange->getBlock("+toString(i)+"))"),
885 *block.range(),("(*block["+toString(i)+","+toString(j)+"].range())")
886 );
887 }
888 if(domainBlock.get()) {
889 THYRA_ASSERT_VEC_SPACES_NAMES(
890 "DefaultBlockedLinearOp<Scalar>::setBlockSpaces(i,j,block):\n\n"
891 "Adding block: " + block.description(),
892 *domainBlock,("(*productDomain->getBlock("+toString(j)+"))"),
893 *block.domain(),("(*block["+toString(i)+","+toString(j)+"].domain())")
894 );
895 }
896#endif // TEUCHOS_DEBUG
897 }
898
899 // Add spaces missing range and domain space blocks if we are doing a
900 // flexible fill (otherwise these loops will not be executed)
901 for( int k = numRowBlocks_; k <= i; ++k )
902 rangeBlocks_.push_back(Teuchos::null);
903 for( int k = numColBlocks_; k <= j; ++k )
904 domainBlocks_.push_back(Teuchos::null);
905
906 // Set the incoming range and domain blocks if not already set
907 if(!productRange_.get()) {
908 if(!rangeBlocks_[i].get())
909 rangeBlocks_[i] = block.range().assert_not_null();
910 if(!domainBlocks_[j].get()) {
911 domainBlocks_[j] = block.domain().assert_not_null();
912 }
913 }
914
915 // Update the current number of row and columns blocks if doing a flexible
916 // fill.
917 if(!Ops_.size()) {
918 numRowBlocks_ = rangeBlocks_.size();
919 numColBlocks_ = domainBlocks_.size();
920 }
921
922}
923
924
925template<class Scalar>
926template<class LinearOpType>
927void DefaultBlockedLinearOp<Scalar>::setBlockImpl(
928 const int i, const int j,
929 const RCP<LinearOpType> &block
930 )
931{
932 setBlockSpaces(i, j, *block);
933 if (Ops_.size()) {
934 // We are doing a fill with a fixed number of row and column blocks so we
935 // can just set this.
936 Ops_[numRowBlocks_*j+i] = block;
937 }
938 else {
939 // We are doing a flexible fill so add the block to the stack of blocks or
940 // replace a block that already exists.
941 bool foundBlock = false;
942 for( unsigned int k = 0; k < Ops_stack_.size(); ++k ) {
943 BlockEntry<Scalar> &block_i_j = Ops_stack_[k];
944 if( block_i_j.i == i && block_i_j.j == j ) {
945 block_i_j.block = block;
946 foundBlock = true;
947 break;
948 }
949 }
950 if(!foundBlock)
951 Ops_stack_.push_back(BlockEntry<Scalar>(i,j,block));
952 }
953}
954
955
956template<class Scalar>
957void DefaultBlockedLinearOp<Scalar>::adjustBlockSpaces()
958{
959
960#ifdef TEUCHOS_DEBUG
961 TEUCHOS_ASSERT_INEQUALITY(Ops_.size(), !=, 0);
962#endif
963
964 //
965 // Loop through the rows and columns looking for rows with mixed
966 // single-space range and/or domain spaces on operators and set the single
967 // spaces as the block space if it exists.
968 //
969 // NOTE: Once we get here, we can safely assume that all of the operators
970 // are compatible w.r.t. their spaces so if there are rows and/or columns
971 // with mixed product and single vector spaces that we can just pick the
972 // single vector space for the whole row and/or column.
973 //
974
975 // Adjust blocks in the range space
976 for (int i = 0; i < numRowBlocks_; ++i) {
977 for (int j = 0; j < numColBlocks_; ++j) {
978 const RCP<const LinearOpBase<Scalar> >
979 op_i_j = Ops_[numRowBlocks_*j+i];
980 if (is_null(op_i_j))
981 continue;
982 const RCP<const VectorSpaceBase<Scalar> > range_i_j = op_i_j->range();
983 if (is_null(productVectorSpaceBase<Scalar>(range_i_j, false))) {
984 rangeBlocks_[i] = range_i_j;
985 break;
986 }
987 }
988 }
989
990 // Adjust blocks in the domain space
991 for (int j = 0; j < numColBlocks_; ++j) {
992 for (int i = 0; i < numRowBlocks_; ++i) {
993 const RCP<const LinearOpBase<Scalar> >
994 op_i_j = Ops_[numRowBlocks_*j+i];
995 if (is_null(op_i_j))
996 continue;
997 const RCP<const VectorSpaceBase<Scalar> >
998 domain_i_j = op_i_j->domain();
999 if (is_null(productVectorSpaceBase<Scalar>(domain_i_j, false))) {
1000 domainBlocks_[j] = domain_i_j;
1001 break;
1002 }
1003 }
1004 }
1005
1006}
1007
1008
1009} // namespace Thyra
1010
1011
1012template<class Scalar>
1013Teuchos::RCP<Thyra::DefaultBlockedLinearOp<Scalar> >
1014Thyra::defaultBlockedLinearOp()
1015{
1016 return Teuchos::rcp(new DefaultBlockedLinearOp<Scalar>());
1017}
1018
1019
1020template<class Scalar>
1021Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1022Thyra::block1x1(
1023 const RCP<const LinearOpBase<Scalar> > &A00,
1024 const std::string &label
1025 )
1026{
1027 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1028 M = defaultBlockedLinearOp<Scalar>();
1029 M->beginBlockFill(1,1);
1030 M->setBlock(0, 0, A00);
1031 M->endBlockFill();
1032 if (label.length())
1033 M->setObjectLabel(label);
1034 return M;
1035}
1036
1037
1038template<class Scalar>
1039Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1040Thyra::block1x2(
1041 const RCP<const LinearOpBase<Scalar> > &A00,
1042 const RCP<const LinearOpBase<Scalar> > &A01,
1043 const std::string &label
1044 )
1045{
1046 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1047 M = defaultBlockedLinearOp<Scalar>();
1048 M->beginBlockFill(1,2);
1049 if(A00.get()) M->setBlock(0,0,A00);
1050 if(A01.get()) M->setBlock(0,1,A01);
1051 M->endBlockFill();
1052 if (label.length())
1053 M->setObjectLabel(label);
1054 return M;
1055}
1056
1057
1058template<class Scalar>
1059Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1060Thyra::block2x1(
1061 const RCP<const LinearOpBase<Scalar> > &A00,
1062 const RCP<const LinearOpBase<Scalar> > &A10,
1063 const std::string &label
1064 )
1065{
1066 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1067 M = defaultBlockedLinearOp<Scalar>();
1068 M->beginBlockFill(2,1);
1069 if(A00.get()) M->setBlock(0,0,A00);
1070 if(A10.get()) M->setBlock(1,0,A10);
1071 M->endBlockFill();
1072 if (label.length())
1073 M->setObjectLabel(label);
1074 return M;
1075}
1076
1077
1078template<class Scalar>
1079Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1080Thyra::block2x2(
1081 const RCP<const LinearOpBase<Scalar> > &A00,
1082 const RCP<const LinearOpBase<Scalar> > &A01,
1083 const RCP<const LinearOpBase<Scalar> > &A10,
1084 const RCP<const LinearOpBase<Scalar> > &A11,
1085 const std::string &label
1086 )
1087{
1088 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1089 M = defaultBlockedLinearOp<Scalar>();
1090 M->beginBlockFill(2,2);
1091 if(A00.get()) M->setBlock(0,0,A00);
1092 if(A01.get()) M->setBlock(0,1,A01);
1093 if(A10.get()) M->setBlock(1,0,A10);
1094 if(A11.get()) M->setBlock(1,1,A11);
1095 M->endBlockFill();
1096 if (label.length())
1097 M->setObjectLabel(label);
1098 return M;
1099}
1100
1101
1102template<class Scalar>
1103Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1104Thyra::nonconstBlock1x1(
1105 const RCP<LinearOpBase<Scalar> > &A00,
1106 const std::string &label
1107 )
1108{
1109 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1110 M = defaultBlockedLinearOp<Scalar>();
1111 M->beginBlockFill(1, 1);
1112 M->setNonconstBlock(0, 0, A00);
1113 M->endBlockFill();
1114 if (label.length())
1115 M->setObjectLabel(label);
1116 return M;
1117}
1118
1119
1120template<class Scalar>
1121Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1122Thyra::nonconstBlock1x2(
1123 const RCP<LinearOpBase<Scalar> > &A00,
1124 const RCP<LinearOpBase<Scalar> > &A01,
1125 const std::string &label
1126 )
1127{
1128 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1129 M = defaultBlockedLinearOp<Scalar>();
1130 M->beginBlockFill(1,2);
1131 if(A00.get()) M->setNonconstBlock(0,0,A00);
1132 if(A01.get()) M->setNonconstBlock(0,1,A01);
1133 M->endBlockFill();
1134 if (label.length())
1135 M->setObjectLabel(label);
1136 return M;
1137}
1138
1139
1140template<class Scalar>
1141Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1142Thyra::nonconstBlock2x1(
1143 const RCP<LinearOpBase<Scalar> > &A00,
1144 const RCP<LinearOpBase<Scalar> > &A10,
1145 const std::string &label
1146 )
1147{
1148 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1149 M = defaultBlockedLinearOp<Scalar>();
1150 M->beginBlockFill(2,1);
1151 if(A00.get()) M->setNonconstBlock(0,0,A00);
1152 if(A10.get()) M->setNonconstBlock(1,0,A10);
1153 M->endBlockFill();
1154 if (label.length())
1155 M->setObjectLabel(label);
1156 return M;
1157}
1158
1159
1160template<class Scalar>
1161Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1162Thyra::nonconstBlock2x2(
1163 const RCP<LinearOpBase<Scalar> > &A00,
1164 const RCP<LinearOpBase<Scalar> > &A01,
1165 const RCP<LinearOpBase<Scalar> > &A10,
1166 const RCP<LinearOpBase<Scalar> > &A11,
1167 const std::string &label
1168 )
1169{
1170 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
1171 M = defaultBlockedLinearOp<Scalar>();
1172 M->beginBlockFill(2,2);
1173 if(A00.get()) M->setNonconstBlock(0,0,A00);
1174 if(A01.get()) M->setNonconstBlock(0,1,A01);
1175 if(A10.get()) M->setNonconstBlock(1,0,A10);
1176 if(A11.get()) M->setNonconstBlock(1,1,A11);
1177 M->endBlockFill();
1178 if (label.length())
1179 M->setObjectLabel(label);
1180 return M;
1181}
1182
1183
1184//
1185// Explicit instantiation macro
1186//
1187// Must be expanded from within the Thyra namespace!
1188//
1189
1190
1191#define THYRA_DEFAULT_BLOCKED_LINEAR_OP_INSTANT(SCALAR) \
1192 \
1193 template class DefaultBlockedLinearOp<SCALAR >; \
1194 \
1195 template RCP<DefaultBlockedLinearOp<SCALAR > > \
1196 defaultBlockedLinearOp<SCALAR >(); \
1197 \
1198 template RCP<const LinearOpBase<SCALAR > > \
1199 block1x1( \
1200 const RCP<const LinearOpBase<SCALAR > > &A00, \
1201 const std::string &label \
1202 ); \
1203 \
1204 template RCP<const LinearOpBase<SCALAR > > \
1205 block1x2( \
1206 const RCP<const LinearOpBase<SCALAR > > &A00, \
1207 const RCP<const LinearOpBase<SCALAR > > &A01, \
1208 const std::string &label \
1209 ); \
1210 \
1211 template RCP<const LinearOpBase<SCALAR > > \
1212 block2x1( \
1213 const RCP<const LinearOpBase<SCALAR > > &A00, \
1214 const RCP<const LinearOpBase<SCALAR > > &A10, \
1215 const std::string &label \
1216 ); \
1217 \
1218 template RCP<const LinearOpBase<SCALAR > > \
1219 block2x2( \
1220 const RCP<const LinearOpBase<SCALAR > > &A00, \
1221 const RCP<const LinearOpBase<SCALAR > > &A01, \
1222 const RCP<const LinearOpBase<SCALAR > > &A10, \
1223 const RCP<const LinearOpBase<SCALAR > > &A11, \
1224 const std::string &label \
1225 ); \
1226 \
1227 template RCP<LinearOpBase<SCALAR > > \
1228 nonconstBlock1x1( \
1229 const RCP<LinearOpBase<SCALAR > > &A00, \
1230 const std::string &label \
1231 ); \
1232 \
1233 template RCP<LinearOpBase<SCALAR > > \
1234 nonconstBlock1x2( \
1235 const RCP<LinearOpBase<SCALAR > > &A00, \
1236 const RCP<LinearOpBase<SCALAR > > &A01, \
1237 const std::string &label \
1238 ); \
1239 \
1240 template RCP<LinearOpBase<SCALAR > > \
1241 nonconstBlock2x1( \
1242 const RCP<LinearOpBase<SCALAR > > &A00, \
1243 const RCP<LinearOpBase<SCALAR > > &A10, \
1244 const std::string &label \
1245 ); \
1246 \
1247 template RCP<LinearOpBase<SCALAR > > \
1248 nonconstBlock2x2( \
1249 const RCP<LinearOpBase<SCALAR > > &A00, \
1250 const RCP<LinearOpBase<SCALAR > > &A01, \
1251 const RCP<LinearOpBase<SCALAR > > &A10, \
1252 const RCP<LinearOpBase<SCALAR > > &A11, \
1253 const std::string &label \
1254 );
1255
1256
1257#endif // THYRA_DEFAULT_BLOCKED_LINEAR_OP_DEF_HPP
Teuchos::Describable::description
virtual std::string description() const
Teuchos::RCP::get
T * get() const
Thyra::DefaultBlockedLinearOp::domain
Teuchos::RCP< const VectorSpaceBase< Scalar > > domain() const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:279
Thyra::DefaultBlockedLinearOp::blockIsConst
bool blockIsConst(const int i, const int j) const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:239
Thyra::DefaultBlockedLinearOp::scaleRightImpl
virtual void scaleRightImpl(const VectorBase< Scalar > &col_scaling)
Definition Thyra_DefaultBlockedLinearOp_def.hpp:740
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:385
Thyra::DefaultBlockedLinearOp::nonconstBlock2x2
Teuchos::RCP< LinearOpBase< Scalar > > nonconstBlock2x2(const Teuchos::RCP< LinearOpBase< Scalar > > &A00, const Teuchos::RCP< LinearOpBase< Scalar > > &A01, const Teuchos::RCP< LinearOpBase< Scalar > > &A10, const Teuchos::RCP< LinearOpBase< Scalar > > &A11, const std::string &label="")
Form an implicit block 2x2 linear operator [ A00, A01; A10, A11 ].
Thyra::DefaultBlockedLinearOp::acceptsBlock
bool acceptsBlock(const int i, const int j) const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:87
Thyra::DefaultBlockedLinearOp::range
Teuchos::RCP< const VectorSpaceBase< Scalar > > range() const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:271
Thyra::DefaultBlockedLinearOp::rowStatIsSupportedImpl
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:460
Thyra::DefaultBlockedLinearOp::getNonconstBlock
Teuchos::RCP< LinearOpBase< Scalar > > getNonconstBlock(const int i, const int j)
Definition Thyra_DefaultBlockedLinearOp_def.hpp:250
Thyra::DefaultBlockedLinearOp::getBlock
Teuchos::RCP< const LinearOpBase< Scalar > > getBlock(const int i, const int j) const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:259
Thyra::DefaultBlockedLinearOp::blockExists
bool blockExists(const int i, const int j) const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:228
Thyra::DefaultBlockedLinearOp::block2x1
Teuchos::RCP< const LinearOpBase< Scalar > > block2x1(const Teuchos::RCP< const LinearOpBase< Scalar > > &A00, const Teuchos::RCP< const LinearOpBase< Scalar > > &A10, const std::string &label="")
Form an implicit block 2x1 linear operator [ A00; A10 ].
Thyra::DefaultBlockedLinearOp::scaleLeftImpl
virtual void scaleLeftImpl(const VectorBase< Scalar > &row_scaling)
Definition Thyra_DefaultBlockedLinearOp_def.hpp:687
Thyra::DefaultBlockedLinearOp::setNonconstBlock
void setNonconstBlock(const int i, const int j, const Teuchos::RCP< LinearOpBase< Scalar > > &block)
Definition Thyra_DefaultBlockedLinearOp_def.hpp:98
Thyra::DefaultBlockedLinearOp::nonconstBlock1x2
Teuchos::RCP< LinearOpBase< Scalar > > nonconstBlock1x2(const Teuchos::RCP< LinearOpBase< Scalar > > &A00, const Teuchos::RCP< LinearOpBase< Scalar > > &A01, const std::string &label="")
Form an implicit block 1x2 linear operator [ A00, A01 ].
Thyra::DefaultBlockedLinearOp::productDomain
Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productDomain() const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:221
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:311
Thyra::DefaultBlockedLinearOp::defaultBlockedLinearOp
RCP< DefaultBlockedLinearOp< Scalar > > defaultBlockedLinearOp()
Nonmember default constructor.
Thyra::DefaultBlockedLinearOp::productRange
Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productRange() const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:213
Thyra::DefaultBlockedLinearOp::getRowStatImpl
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Teuchos::Ptr< VectorBase< Scalar > > &rowStatVec) const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:525
Thyra::DefaultBlockedLinearOp::setBlock
void setBlock(const int i, const int j, const Teuchos::RCP< const LinearOpBase< Scalar > > &block)
Definition Thyra_DefaultBlockedLinearOp_def.hpp:108
Thyra::DefaultBlockedLinearOp::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:297
Thyra::DefaultBlockedLinearOp::block1x2
Teuchos::RCP< const LinearOpBase< Scalar > > block1x2(const Teuchos::RCP< const LinearOpBase< Scalar > > &A00, const Teuchos::RCP< const LinearOpBase< Scalar > > &A01, const std::string &label="")
Form an implicit block 1x2 linear operator [ A00, A01 ].
Thyra::DefaultBlockedLinearOp::nonconstBlock2x1
Teuchos::RCP< LinearOpBase< Scalar > > nonconstBlock2x1(const Teuchos::RCP< LinearOpBase< Scalar > > &A00, const Teuchos::RCP< LinearOpBase< Scalar > > &A10, const std::string &label="")
Form an implicit block 2x1 linear operator [ A00; A10 ].
Thyra::DefaultBlockedLinearOp::block2x2
Teuchos::RCP< const LinearOpBase< Scalar > > block2x2(const Teuchos::RCP< const LinearOpBase< Scalar > > &A00, const Teuchos::RCP< const LinearOpBase< Scalar > > &A01, const Teuchos::RCP< const LinearOpBase< Scalar > > &A10, const Teuchos::RCP< const LinearOpBase< Scalar > > &A11, const std::string &label="")
Form an implicit block 2x2 linear operator [ A00, A01; A10, A11 ].
Thyra::DefaultBlockedLinearOp::clone
Teuchos::RCP< const LinearOpBase< Scalar > > clone() const
Definition Thyra_DefaultBlockedLinearOp_def.hpp:287
Thyra::DefaultBlockedLinearOp::opSupportedImpl
bool opSupportedImpl(EOpTransp M_trans) const
Returns true only if all constituent operators support M_trans.
Definition Thyra_DefaultBlockedLinearOp_def.hpp:369
Thyra::DefaultBlockedLinearOp::block1x1
Teuchos::RCP< const LinearOpBase< Scalar > > block1x1(const Teuchos::RCP< const LinearOpBase< Scalar > > &A00, const std::string &label="")
Form an implicit block 1x1 linear operator [ A00 ].
Thyra::DefaultBlockedLinearOp::nonconstBlock1x1
Teuchos::RCP< LinearOpBase< Scalar > > nonconstBlock1x1(const Teuchos::RCP< LinearOpBase< Scalar > > &A00, const std::string &label="")
Form an implicit block 1x1 linear operator [ A00 ].
Thyra::DefaultProductMultiVector::castOrCreateSingleBlockProductMultiVector
RCP< const ProductMultiVectorBase< Scalar > > castOrCreateSingleBlockProductMultiVector(const RCP< const DefaultProductVectorSpace< Scalar > > &productSpace, const RCP< const MultiVectorBase< Scalar > > &mv)
Dynamic cast to a const product multi-vector or create a new product multi-vector with one component ...
Thyra::DefaultProductMultiVector::nonconstCastOrCreateSingleBlockProductMultiVector
RCP< ProductMultiVectorBase< Scalar > > nonconstCastOrCreateSingleBlockProductMultiVector(const RCP< const DefaultProductVectorSpace< Scalar > > &productSpace, const RCP< MultiVectorBase< Scalar > > &mv)
Dynamic cast to a const product multi-vector or create a new product multi-vector with one component ...
Thyra::DefaultProductVectorSpace::productVectorSpace
RCP< DefaultProductVectorSpace< Scalar > > productVectorSpace()
Nonmember constructor that constructs to uninitialized.
Definition Thyra_DefaultProductVectorSpace_decl.hpp:361
Thyra::DefaultScaledAdjointLinearOp::scale
RCP< const LinearOpBase< Scalar > > scale(const Scalar &scalar, const RCP< const LinearOpBase< Scalar > > &Op, const std::string &label="")
Build an implicit const scaled linear operator.
Thyra::LinearOpBase
Base class for all linear operators.
Definition Thyra_LinearOpBase_decl.hpp:159
Thyra::LinearOpBase::apply
void apply(const LinearOpBase< Scalar > &M, const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha=static_cast< Scalar >(1.0), const Scalar beta=static_cast< Scalar >(0.0))
Non-member function call for M.apply(...).
Definition Thyra_LinearOpBase_def.hpp:41
Thyra::LinearOpBase::opSupported
bool opSupported(const LinearOpBase< Scalar > &M, EOpTransp M_trans)
Determines if an operation is supported for a single scalar type.
Definition Thyra_LinearOpBase_decl.hpp:431
Thyra::MultiVectorBase
Interface for a collection of column vectors called a multi-vector.
Definition Thyra_MultiVectorBase_decl.hpp:464
Thyra::MultiVectorBase::Vp_V
void Vp_V(const Ptr< MultiVectorBase< Scalar > > &Z, const MultiVectorBase< Scalar > &X)
Z(i,j) += X(i,j), i = 0...Z->range()->dim()-1, j = 0...Z->domain()->dim()-1.
Definition Thyra_MultiVectorStdOps_def.hpp:228
Thyra::ProductVectorBase
Base interface for product vectors.
Definition Thyra_ProductVectorBase.hpp:58
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::ProductVectorSpaceBase::productVectorSpaceBase
RCP< const ProductVectorSpaceBase< Scalar > > productVectorSpaceBase(const RCP< const VectorSpaceBase< Scalar > > &v, const bool forceSuccess=true)
Dynamic cast from a const VectorSpaceBase to a const ProductVectorSpaceBase object and thow exception...
Definition Thyra_ProductVectorSpaceBase.hpp:122
Thyra::RowStatLinearOpBase
Interface for exxtracting row statistics as a VectorBase from a supporting LinearOpBase object.
Definition Thyra_RowStatLinearOpBase.hpp:44
Thyra::ScaledLinearOpBase
Applies left or right sclaing to the linear operator.
Definition Thyra_ScaledLinearOpBase.hpp:50
Thyra::VectorBase
Abstract interface for finite-dimensional dense vectors.
Definition Thyra_VectorBase.hpp:115
Thyra::VectorBase::put_scalar
void put_scalar(const Scalar &alpha, const Ptr< VectorBase< Scalar > > &y)
Assign all elements to a scalar: y(i) = alpha, i = 0...y->space()->dim()-1.
Definition Thyra_VectorStdOps_def.hpp:135
Thyra::VectorBase::reciprocal
void reciprocal(const VectorBase< Scalar > &x, const Ptr< VectorBase< Scalar > > &y)
Element-wise reciprocal: y(i) = 1/x(i), i = 0...y->space()->dim()-1.
Definition Thyra_VectorStdOps_def.hpp:175
Thyra::VectorSpaceBase::createMember
RCP< VectorBase< Scalar > > createMember(const RCP< const VectorSpaceBase< Scalar > > &vs, const std::string &label="")
Create a vector member from the vector space.
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_ASSERT_INEQUALITY
#define TEUCHOS_ASSERT_INEQUALITY(val1, comp, val2)
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES
#define THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the multi-vector vers...
Definition Thyra_AssertOp.hpp:218
THYRA_ASSERT_VEC_SPACES_NAMES
#define THYRA_ASSERT_VEC_SPACES_NAMES(FUNC_NAME, VS1, VS1_NAME, VS2, VS2_NAME)
Helper assertion macro.
Definition Thyra_AssertOp.hpp:134
Thyra::unwrap
void unwrap(const LinearOpBase< Scalar > &Op, Scalar *scalar, EOpTransp *transp, const LinearOpBase< Scalar > **origOp)
Extract the overallScalar, overallTransp and const origOp from a const LinearOpBase object.
Definition Thyra_ScaledAdjointLinearOpBase_def.hpp:18
Thyra::toString
const char * toString(EOpTransp transp)
Return a string name for a EOpTransp value.
Definition Thyra_OperatorVectorTypes.hpp:149
Thyra::NOTRANS
NOTRANS
Type for the dimension of a vector space. `**/ typedef Teuchos::Ordinal Ordinal;.
Definition Thyra_OperatorVectorTypes.hpp:130
Thyra::real_trans
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
Definition Thyra_OperatorVectorTypes.hpp:168
Teuchos::as
TypeTo as(const TypeFrom &t)
Teuchos::dyn_cast
T_To & dyn_cast(T_From &from)
Teuchos::toString
std::string toString(const T &t)
Teuchos::FancyOStream
basic_FancyOStream< char > FancyOStream
Teuchos::OSTab
basic_OSTab< char > OSTab
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::VERB_MEDIUM
VERB_MEDIUM
Teuchos::VERB_HIGH
VERB_HIGH
Teuchos::VERB_EXTREME
VERB_EXTREME
Teuchos::VERB_DEFAULT
VERB_DEFAULT
Teuchos::VERB_LOW
VERB_LOW
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Generated by doxygen 1.15.0