42 #ifndef THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP 43 #define THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP 46 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolve_decl.hpp" 47 #include "Thyra_ProductMultiVectorBase.hpp" 48 #include "Thyra_DefaultProductVectorSpace.hpp" 49 #include "Thyra_AssertOp.hpp" 61 template<
class Scalar>
63 : blockFillIsActive_(false), numDiagBlocks_(0)
67 template<
class Scalar>
72 assertAndSetBlockStructure(*blocks);
73 blocks_.initialize(blocks);
77 template<
class Scalar>
82 assertAndSetBlockStructure(*blocks);
83 blocks_.initialize(blocks);
87 template<
class Scalar>
88 RCP<PhysicallyBlockedLinearOpBase<Scalar> >
91 return blocks_.getNonconstObj();
95 template<
class Scalar>
96 RCP<const PhysicallyBlockedLinearOpBase<Scalar> >
99 return blocks_.getConstObj();
106 template<
class Scalar>
108 const int i,
const int j
111 assertBlockFillIsActive(
true);
112 assertBlockRowCol(i,j);
116 template<
class Scalar>
118 const int i,
const int j,
122 setLOWSBlockImpl(i,j,block);
126 template<
class Scalar>
128 const int i,
const int j,
132 setLOWSBlockImpl(i,j,block);
139 template<
class Scalar>
142 assertBlockFillIsActive(
false);
143 TEUCHOS_TEST_FOR_EXCEPT(
"ToDo: Have not implemented flexible block fill yet!");
147 template<
class Scalar>
149 const int numRowBlocks,
const int numColBlocks
154 TEUCHOS_ASSERT_EQUALITY(numRowBlocks, numColBlocks);
156 assertBlockFillIsActive(
false);
157 numDiagBlocks_ = numRowBlocks;
158 diagonalBlocks_.resize(numDiagBlocks_);
159 productRange_ = null;
160 productDomain_ = null;
161 blockFillIsActive_ =
true;
165 template<
class Scalar>
172 TEUCHOS_TEST_FOR_EXCEPT( is_null(productRange_in) );
173 TEUCHOS_TEST_FOR_EXCEPT( is_null(productDomain_in) );
174 TEUCHOS_TEST_FOR_EXCEPT( productRange_in->numBlocks() != productDomain_in->numBlocks() );
176 assertBlockFillIsActive(
false);
177 productRange_ = productRange_in;
178 productDomain_ = productDomain_in;
179 numDiagBlocks_ = productRange_in->numBlocks();
180 diagonalBlocks_.resize(numDiagBlocks_);
181 blockFillIsActive_ =
true;
185 template<
class Scalar>
188 return blockFillIsActive_;
192 template<
class Scalar>
194 const int i,
const int j
197 assertBlockFillIsActive(
true);
198 assertBlockRowCol(i,j);
203 template<
class Scalar>
205 const int i,
const int j,
209 assertBlockFillIsActive(
true);
210 TEUCHOS_TEST_FOR_EXCEPT(
"Error, we don't support off-diagonal LOB objects yet!");
214 template<
class Scalar>
216 const int i,
const int j,
220 assertBlockFillIsActive(
true);
221 TEUCHOS_TEST_FOR_EXCEPT(
"Error, we don't support off-diagonal LOB objects yet!");
225 template<
class Scalar>
228 assertBlockFillIsActive(
true);
229 Array<RCP<const VectorSpaceBase<Scalar> > > rangeSpaces;
230 Array<RCP<const VectorSpaceBase<Scalar> > > domainSpaces;
231 for (
int k = 0; k < numDiagBlocks_; ++k ) {
232 const RCP<const LinearOpWithSolveBase<Scalar> > lows_k =
233 diagonalBlocks_[k].getConstObj();
234 TEUCHOS_TEST_FOR_EXCEPTION(is_null(lows_k), std::logic_error,
235 "Error, the block diagonal k="<<k<<
" can not be null when ending block fill!" 237 if (is_null(productRange_)) {
238 rangeSpaces.push_back(lows_k->range());
239 domainSpaces.push_back(lows_k->domain());
242 if (is_null(productRange_)) {
243 productRange_ = productVectorSpace<Scalar>(rangeSpaces());
244 productDomain_ = productVectorSpace<Scalar>(domainSpaces());
246 blockFillIsActive_ =
false;
250 template<
class Scalar>
253 assertBlockFillIsActive(
false);
254 productRange_ = Teuchos::null;
255 productDomain_ = Teuchos::null;
257 diagonalBlocks_.resize(0);
264 template<
class Scalar>
265 Teuchos::RCP<LinearOpWithSolveBase<Scalar> >
267 const int i,
const int j
270 assertBlockFillIsActive(
false);
271 assertBlockRowCol(i,j);
273 return Teuchos::null;
274 return diagonalBlocks_[i].getNonconstObj();
278 template<
class Scalar>
279 Teuchos::RCP<const LinearOpWithSolveBase<Scalar> >
281 const int i,
const int j
284 assertBlockFillIsActive(
false);
285 assertBlockRowCol(i, j);
287 return Teuchos::null;
288 return diagonalBlocks_[i].getConstObj();
295 template<
class Scalar>
296 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
299 return productRange_;
303 template<
class Scalar>
304 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
307 return productDomain_;
311 template<
class Scalar>
313 const int i,
const int j
316 assertBlockFillIsActive(
false);
317 assertBlockRowCol(i,j);
320 return !is_null(diagonalBlocks_[i].getConstObj());
324 template<
class Scalar>
326 const int i,
const int j
329 assertBlockFillIsActive(
true);
330 assertBlockRowCol(i,j);
331 return diagonalBlocks_[i].isConst();
335 template<
class Scalar>
336 Teuchos::RCP<LinearOpBase<Scalar> >
338 const int i,
const int j
341 assertBlockFillIsActive(
true);
342 assertBlockRowCol(i,j);
344 return Teuchos::null;
345 return this->getNonconstLOWSBlock(i,j);
349 template<
class Scalar>
350 Teuchos::RCP<const LinearOpBase<Scalar> >
352 const int i,
const int j
355 assertBlockFillIsActive(
true);
356 assertBlockRowCol(i,j);
358 return Teuchos::null;
359 return this->getLOWSBlock(i,j);
366 template<
class Scalar>
367 Teuchos::RCP< const VectorSpaceBase<Scalar> >
370 return this->productRange();
374 template<
class Scalar>
375 Teuchos::RCP< const VectorSpaceBase<Scalar> >
378 return this->productDomain();
382 template<
class Scalar>
383 Teuchos::RCP<const LinearOpBase<Scalar> >
386 return Teuchos::null;
393 template<
class Scalar>
397 assertBlockFillIsActive(
false);
398 std::ostringstream oss;
400 << Teuchos::Describable::description() <<
"{" 401 <<
"numDiagBlocks="<<numDiagBlocks_
407 template<
class Scalar>
409 Teuchos::FancyOStream &out,
410 const Teuchos::EVerbosityLevel verbLevel
413 assertBlockFillIsActive(
false);
414 Teuchos::Describable::describe(out, verbLevel);
425 template<
class Scalar>
430 using Thyra::opSupported;
431 assertBlockFillIsActive(
false);
432 for (
int k = 0; k < numDiagBlocks_; ++k ) {
433 if ( !opSupported(*diagonalBlocks_[k].getConstObj(),M_trans) )
443 template<
class Scalar>
454 using Teuchos::dyn_cast;
459 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
460 *
this, M_trans, X_in, &*Y_inout
462 #endif // THYRA_DEBUG 479 for (
int i = 0; i < numDiagBlocks_; ++ i ) {
480 Thyra::apply( *diagonalBlocks_[i].getConstObj(), M_trans,
492 template<
class Scalar>
498 assertBlockFillIsActive(
false);
499 for (
int k = 0; k < numDiagBlocks_; ++k ) {
500 if ( !Thyra::solveSupports( *diagonalBlocks_[k].getConstObj(), M_trans ) )
510 template<
class Scalar>
516 using Thyra::solveSupportsSolveMeasureType;
517 assertBlockFillIsActive(
false);
518 for (
int k = 0; k < numDiagBlocks_; ++k ) {
520 !solveSupportsSolveMeasureType(
521 *diagonalBlocks_[k].getConstObj(),
522 M_trans, solveMeasureType
533 template<
class Scalar>
544 using Teuchos::dyn_cast;
549 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
550 *
this, M_trans, *X_inout, &B_in
552 TEUCHOS_TEST_FOR_EXCEPT(!this->solveSupportsImpl(M_trans));
553 TEUCHOS_TEST_FOR_EXCEPTION(
554 nonnull(solveCriteria) && !solveCriteria->solveMeasureType.useDefault(),
556 "Error, we can't handle any non-default solve criteria yet!" 560 #endif // THYRA_DEBUG 577 for (
int i = 0; i < numDiagBlocks_; ++ i ) {
578 const RCP<const LinearOpWithSolveBase<Scalar> >
579 Op_k = diagonalBlocks_[i].getConstObj();
580 Op_k->setOStream(this->getOStream());
581 Op_k->setVerbLevel(this->getVerbLevel());
582 Thyra::solve( *Op_k, M_trans, *B.getMultiVectorBlock(i),
583 X.getNonconstMultiVectorBlock(i).ptr() );
596 template<
class Scalar>
598 bool blockFillIsActive_in
602 TEUCHOS_TEST_FOR_EXCEPT(!(blockFillIsActive_==blockFillIsActive_in));
607 template<
class Scalar>
608 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertBlockRowCol(
609 const int i,
const int j
613 TEUCHOS_TEST_FOR_EXCEPTION(
614 !( 0 <= i && i < numDiagBlocks_ ), std::logic_error,
615 "Error, i="<<i<<
" does not fall in the range [0,"<<numDiagBlocks_-1<<
"]!" 617 TEUCHOS_TEST_FOR_EXCEPTION(
618 !( 0 <= j && j < numDiagBlocks_ ), std::logic_error,
619 "Error, j="<<j<<
" does not fall in the range [0,"<<numDiagBlocks_-1<<
"]!" 625 template<
class Scalar>
626 template<
class LinearOpWithSolveType>
627 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setLOWSBlockImpl(
628 const int i,
const int j,
629 const Teuchos::RCP<LinearOpWithSolveType> &block
632 assertBlockFillIsActive(
true);
634 TEUCHOS_ASSERT_INEQUALITY( i, >=, 0 );
635 TEUCHOS_ASSERT_INEQUALITY( j, >=, 0 );
636 TEUCHOS_ASSERT_INEQUALITY( i, <, numDiagBlocks_ );
637 TEUCHOS_ASSERT_INEQUALITY( j, <, numDiagBlocks_ );
638 TEUCHOS_TEST_FOR_EXCEPTION(
639 !this->acceptsLOWSBlock(i,j), std::logic_error,
640 "Error, this DefaultBlockedTriangularLinearOpWithSolve does not accept\n" 641 "LOWSB objects for the block i="<<i<<
", j="<<j<<
"!" 644 diagonalBlocks_[i] = block;
648 template<
class Scalar>
649 void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(
650 const PhysicallyBlockedLinearOpBase<Scalar>& blocks
655 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
656 *blocks.range(), *this->range()
659 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
660 *blocks.domain(), *this->domain()
672 #endif // THYRA_DEFAULT_BLOCKED_TRIANGULAR_LINEAR_OP_WITH_SOLVE_DEF_HPP Base interface for product multi-vectors.
RCP< LinearOpWithSolveBase< Scalar > > getNonconstLOWSBlock(const int i, const int j)
#define THYRA_ASSERT_VEC_SPACES(FUNC_NAME, VS1, VS2)
This is a very useful macro that should be used to validate that two vector spaces are compatible...
void setBlocks(const RCP< const PhysicallyBlockedLinearOpBase< Scalar > > &blocks)
Base class for all linear operators that can support a high-level solve operation.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
void setNonconstBlock(const int i, const int j, const RCP< LinearOpBase< Scalar > > &block)
Concrete composite LinearOpWithSolveBase subclass that creates single upper or lower block triangular...
#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...
RCP< const LinearOpBase< Scalar > > clone() const
bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
void setBlock(const int i, const int j, const RCP< const LinearOpBase< Scalar > > &block)
RCP< const PhysicallyBlockedLinearOpBase< Scalar > > getBlocks()
RCP< const VectorSpaceBase< Scalar > > domain() const
SolveStatus< Scalar > solveImpl(const EOpTransp transp, const MultiVectorBase< Scalar > &B, const Ptr< MultiVectorBase< Scalar > > &X, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
bool blockFillIsActive() const
RCP< const LinearOpWithSolveBase< Scalar > > getLOWSBlock(const int i, const int j) const
virtual Teuchos::RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const =0
Returns a non-persisting const view of the (zero-based) kth block multi-vector.
bool blockExists(const int i, const int j) const
RCP< const LinearOpBase< Scalar > > getBlock(const int i, const int j) const
RCP< LinearOpBase< Scalar > > getNonconstBlock(const int i, const int j)
RCP< const ProductVectorSpaceBase< Scalar > > productRange() const
bool opSupportedImpl(EOpTransp M_trans) const
Interface for a collection of column vectors called a multi-vector.
RCP< const VectorSpaceBase< Scalar > > range() const
void setLOWSBlock(const int i, const int j, const RCP< const LinearOpWithSolveBase< Scalar > > &block)
RCP< PhysicallyBlockedLinearOpBase< Scalar > > getNonconstBlocks()
Base interface for physically blocked linear operators.
RCP< const ProductVectorSpaceBase< Scalar > > productDomain() const
Simple struct for the return status from a solve.
Base class for all linear operators.
std::string description() const
Prints just the name DefaultBlockedTriangularLinearOpWithSolve along with the overall dimensions and ...
void setNonconstLOWSBlock(const int i, const int j, const RCP< LinearOpWithSolveBase< Scalar > > &block)
bool blockIsConst(const int i, const int j) const
bool solveSupportsImpl(EOpTransp M_trans) const
bool acceptsBlock(const int i, const int j) const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent linear operators.
Simple struct that defines the requested solution criteria for a solve.
bool acceptsLOWSBlock(const int i, const int j) const
DefaultBlockedTriangularLinearOpWithSolve()
void setNonconstBlocks(const RCP< PhysicallyBlockedLinearOpBase< Scalar > > &blocks)