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"
63 : blockFillIsActive_(false), numDiagBlocks_(0)
72 assertAndSetBlockStructure(*blocks);
73 blocks_.initialize(blocks);
82 assertAndSetBlockStructure(*blocks);
83 blocks_.initialize(blocks);
91 return blocks_.getNonconstObj();
99 return blocks_.getConstObj();
106template<
class Scalar>
108 const int i,
const int j
111 assertBlockFillIsActive(
true);
112 assertBlockRowCol(i,j);
116template<
class Scalar>
118 const int i,
const int j,
122 setLOWSBlockImpl(i,j,block);
126template<
class Scalar>
128 const int i,
const int j,
132 setLOWSBlockImpl(i,j,block);
139template<
class Scalar>
142 assertBlockFillIsActive(
false);
147template<
class Scalar>
149 const int numRowBlocks,
const int numColBlocks
158 assertBlockFillIsActive(
false);
159 numDiagBlocks_ = numRowBlocks;
160 diagonalBlocks_.resize(numDiagBlocks_);
161 productRange_ = null;
162 productDomain_ = null;
163 blockFillIsActive_ =
true;
167template<
class Scalar>
178 assertBlockFillIsActive(
false);
179 productRange_ = productRange_in;
180 productDomain_ = productDomain_in;
181 numDiagBlocks_ = productRange_in->numBlocks();
182 diagonalBlocks_.resize(numDiagBlocks_);
183 blockFillIsActive_ =
true;
187template<
class Scalar>
190 return blockFillIsActive_;
194template<
class Scalar>
196 const int i,
const int j
199 assertBlockFillIsActive(
true);
200 assertBlockRowCol(i,j);
205template<
class Scalar>
207 const int ,
const int ,
211 assertBlockFillIsActive(
true);
216template<
class Scalar>
218 const int ,
const int ,
222 assertBlockFillIsActive(
true);
227template<
class Scalar>
230 assertBlockFillIsActive(
true);
233 for (
int k = 0; k < numDiagBlocks_; ++k ) {
235 diagonalBlocks_[k].getConstObj();
237 "Error, the block diagonal k="<<k<<
" can not be null when ending block fill!"
239 if (is_null(productRange_)) {
241 domainSpaces.
push_back(lows_k->domain());
244 if (is_null(productRange_)) {
245 productRange_ = productVectorSpace<Scalar>(rangeSpaces());
246 productDomain_ = productVectorSpace<Scalar>(domainSpaces());
248 blockFillIsActive_ =
false;
252template<
class Scalar>
255 assertBlockFillIsActive(
false);
256 productRange_ = Teuchos::null;
257 productDomain_ = Teuchos::null;
259 diagonalBlocks_.resize(0);
266template<
class Scalar>
269 const int i,
const int j
272 assertBlockFillIsActive(
false);
273 assertBlockRowCol(i,j);
275 return Teuchos::null;
276 return diagonalBlocks_[i].getNonconstObj();
280template<
class Scalar>
283 const int i,
const int j
286 assertBlockFillIsActive(
false);
287 assertBlockRowCol(i, j);
289 return Teuchos::null;
290 return diagonalBlocks_[i].getConstObj();
297template<
class Scalar>
301 return productRange_;
305template<
class Scalar>
309 return productDomain_;
313template<
class Scalar>
315 const int i,
const int j
318 assertBlockFillIsActive(
false);
319 assertBlockRowCol(i,j);
322 return !is_null(diagonalBlocks_[i].getConstObj());
326template<
class Scalar>
328 const int i,
const int j
331 assertBlockFillIsActive(
true);
332 assertBlockRowCol(i,j);
333 return diagonalBlocks_[i].isConst();
337template<
class Scalar>
340 const int i,
const int j
343 assertBlockFillIsActive(
true);
344 assertBlockRowCol(i,j);
346 return Teuchos::null;
347 return this->getNonconstLOWSBlock(i,j);
351template<
class Scalar>
354 const int i,
const int j
357 assertBlockFillIsActive(
true);
358 assertBlockRowCol(i,j);
360 return Teuchos::null;
361 return this->getLOWSBlock(i,j);
368template<
class Scalar>
372 return this->productRange();
376template<
class Scalar>
380 return this->productDomain();
384template<
class Scalar>
388 return Teuchos::null;
395template<
class Scalar>
399 assertBlockFillIsActive(
false);
400 std::ostringstream oss;
403 <<
"numDiagBlocks="<<numDiagBlocks_
409template<
class Scalar>
415 assertBlockFillIsActive(
false);
427template<
class Scalar>
432 using Thyra::opSupported;
433 assertBlockFillIsActive(
false);
434 for (
int k = 0; k < numDiagBlocks_; ++k ) {
435 if ( !opSupported(*diagonalBlocks_[k].getConstObj(),M_trans) )
445template<
class Scalar>
461 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
462 *
this, M_trans, X_in, &*Y_inout
477 &X = dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in);
479 &Y = dyn_cast<ProductMultiVectorBase<Scalar> >(*Y_inout);
481 for (
int i = 0; i < numDiagBlocks_; ++ i ) {
482 Thyra::apply( *diagonalBlocks_[i].getConstObj(), M_trans,
494template<
class Scalar>
500 assertBlockFillIsActive(
false);
501 for (
int k = 0; k < numDiagBlocks_; ++k ) {
502 if ( !Thyra::solveSupports( *diagonalBlocks_[k].getConstObj(), M_trans ) )
512template<
class Scalar>
518 using Thyra::solveSupportsSolveMeasureType;
519 assertBlockFillIsActive(
false);
520 for (
int k = 0; k < numDiagBlocks_; ++k ) {
522 !solveSupportsSolveMeasureType(
523 *diagonalBlocks_[k].getConstObj(),
524 M_trans, solveMeasureType
535template<
class Scalar>
551 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::apply(...)",
552 *
this, M_trans, *X_inout, &B_in
575 &B = dyn_cast<const ProductMultiVectorBase<Scalar> >(B_in);
577 &X = dyn_cast<ProductMultiVectorBase<Scalar> >(*X_inout);
579 bool converged =
true;
580 for (
int i = 0; i < numDiagBlocks_; ++ i ) {
582 Op_k = diagonalBlocks_[i].getConstObj();
583 Op_k->setOStream(this->getOStream());
584 Op_k->setVerbLevel(this->getVerbLevel());
587 X.getNonconstMultiVectorBlock(i).ptr(), solveCriteria );
606template<
class Scalar>
608 bool blockFillIsActive_in
614 (void)blockFillIsActive_in;
619template<
class Scalar>
620void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertBlockRowCol(
621 const int i,
const int j
626 !( 0 <= i && i < numDiagBlocks_ ), std::logic_error,
627 "Error, i="<<i<<
" does not fall in the range [0,"<<numDiagBlocks_-1<<
"]!"
630 !( 0 <= j && j < numDiagBlocks_ ), std::logic_error,
631 "Error, j="<<j<<
" does not fall in the range [0,"<<numDiagBlocks_-1<<
"]!"
640template<
class Scalar>
641template<
class LinearOpWithSolveType>
642void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::setLOWSBlockImpl(
643 const int i,
const int j,
647 assertBlockFillIsActive(
true);
654 !this->acceptsLOWSBlock(i,j), std::logic_error,
655 "Error, this DefaultBlockedTriangularLinearOpWithSolve does not accept\n"
656 "LOWSB objects for the block i="<<i<<
", j="<<j<<
"!"
661 diagonalBlocks_[i] = block;
665template<
class Scalar>
666void DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(
667 const PhysicallyBlockedLinearOpBase<Scalar>& blocks
672 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
673 *blocks.range(), *this->range()
676 "DefaultBlockedTriangularLinearOpWithSolve<Scalar>::assertAndSetBlockStructure(blocks)",
677 *blocks.domain(), *this->domain()
void push_back(const value_type &x)
virtual void describe(FancyOStream &out, const EVerbosityLevel verbLevel=verbLevel_default) const
virtual std::string description() const
Concrete composite LinearOpWithSolveBase subclass that creates single upper or lower block triangular...
bool acceptsLOWSBlock(const int i, const int j) const
RCP< const PhysicallyBlockedLinearOpBase< Scalar > > getBlocks()
RCP< const ProductVectorSpaceBase< Scalar > > productDomain() const
void setNonconstBlocks(const RCP< PhysicallyBlockedLinearOpBase< Scalar > > &blocks)
SolveStatus< Scalar > solveImpl(const EOpTransp transp, const MultiVectorBase< Scalar > &B, const Ptr< MultiVectorBase< Scalar > > &X, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
void setLOWSBlock(const int i, const int j, const RCP< const LinearOpWithSolveBase< Scalar > > &block)
RCP< const LinearOpWithSolveBase< Scalar > > getLOWSBlock(const int i, const int j) const
RCP< const LinearOpBase< Scalar > > getBlock(const int i, const int j) const
bool opSupportedImpl(EOpTransp M_trans) const
void setNonconstBlock(const int i, const int j, const RCP< LinearOpBase< Scalar > > &block)
bool solveSupportsImpl(EOpTransp M_trans) const
void setBlocks(const RCP< const PhysicallyBlockedLinearOpBase< Scalar > > &blocks)
bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
RCP< PhysicallyBlockedLinearOpBase< Scalar > > getNonconstBlocks()
RCP< const VectorSpaceBase< Scalar > > domain() const
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)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent linear operators.
void setBlock(const int i, const int j, const RCP< const LinearOpBase< Scalar > > &block)
DefaultBlockedTriangularLinearOpWithSolve()
bool blockIsConst(const int i, const int j) const
RCP< const VectorSpaceBase< Scalar > > range() const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
RCP< const LinearOpBase< Scalar > > clone() const
bool acceptsBlock(const int i, const int j) const
RCP< LinearOpWithSolveBase< Scalar > > getNonconstLOWSBlock(const int i, const int j)
RCP< LinearOpBase< Scalar > > getNonconstBlock(const int i, const int j)
bool blockFillIsActive() const
RCP< const ProductVectorSpaceBase< Scalar > > productRange() const
bool blockExists(const int i, const int j) const
Base class for all linear operators.
Base class for all linear operators that can support a high-level solve operation.
Interface for a collection of column vectors called a multi-vector.
Base interface for physically blocked linear operators.
Base interface for product multi-vectors.
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.
#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)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
@ SOLVE_STATUS_UNCONVERGED
The requested solution criteria has likely not been achieved.
@ SOLVE_STATUS_CONVERGED
The requested solution criteria has likely been achieved.
#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.
#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...
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
T_To & dyn_cast(T_From &from)
Simple struct that defines the requested solution criteria for a solve.
Simple struct for the return status from a solve.
ESolveStatus solveStatus
The return status of the solve.