47 #include "Teko_BlockLowerTriInverseOp.hpp"
49 #include "Teuchos_Utils.hpp"
55 BlockLowerTriInverseOp::BlockLowerTriInverseOp(BlockedLinearOp& L,
56 const std::vector<LinearOp>& invDiag)
61 int blocks = blockRowCount(
L_);
62 TEUCHOS_ASSERT(blocks > 0);
63 TEUCHOS_ASSERT(blocks == blockColCount(
L_));
64 TEUCHOS_ASSERT(blocks == (
int)
invDiag_.size());
87 const double alpha,
const double beta)
const {
88 int blocks = blockCount(src);
90 TEUCHOS_ASSERT(blocks == blockRowCount(
L_));
91 TEUCHOS_ASSERT(blocks == blockCount(dst));
94 srcScrap_ = datacopy(src, srcScrap_);
95 BlockedMultiVector dstCopy;
97 dstScrap_ = datacopy(dst, dstScrap_);
104 std::vector<MultiVector> dstVec;
105 std::vector<MultiVector> scrapVec;
106 for (
int b = 0; b < blocks; b++) {
107 dstVec.push_back(getBlock(b, dstCopy));
108 scrapVec.push_back(getBlock(b, srcScrap_));
113 for (
int b = 0; b < blocks; b++) {
114 applyOp(
invDiag_[b], scrapVec[b], dstVec[b]);
117 for (
int i = b + 1; i < blocks; i++) {
118 LinearOp u_ib = getBlock(i, b,
L_);
119 if (u_ib != Teuchos::null) {
120 applyOp(u_ib, dstVec[b], scrapVec[i], -1.0, 1.0);
127 update(alpha, dstCopy, beta, dst);
128 else if (alpha != 1.0)
132 void BlockLowerTriInverseOp::describe(Teuchos::FancyOStream& out_arg,
133 const Teuchos::EVerbosityLevel verbLevel)
const {
134 using Teuchos::OSTab;
136 RCP<Teuchos::FancyOStream> out = rcp(&out_arg,
false);
139 case Teuchos::VERB_DEFAULT:
140 case Teuchos::VERB_LOW: *out << this->description() << std::endl;
break;
141 case Teuchos::VERB_MEDIUM:
142 case Teuchos::VERB_HIGH:
143 case Teuchos::VERB_EXTREME: {
144 *out << Teuchos::Describable::description() <<
"{"
145 <<
"rangeDim=" << this->
range()->dim() <<
",domainDim=" << this->
domain()->dim()
146 <<
",rows=" << blockRowCount(
L_) <<
",cols=" << blockColCount(
L_) <<
"}\n";
149 *out <<
"[L Operator] = ";
150 *out << Teuchos::describe(*
L_, verbLevel);
154 *out <<
"[invDiag Operators]:\n";
156 for (
int i = 0; i < blockRowCount(
L_); i++) {
157 *out <<
"[invD(" << i <<
")] = ";
158 *out << Teuchos::describe(*
invDiag_[i], verbLevel);
163 default: TEUCHOS_TEST_FOR_EXCEPT(
true);
virtual VectorSpace domain() const
Domain space of this operator.
const BlockedLinearOp L_
operator
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productDomain_
Domain vector space.
virtual void implicitApply(const BlockedMultiVector &x, BlockedMultiVector &y, const double alpha=1.0, const double beta=0.0) const
Perform a matrix vector multiply with this operator.
std::vector< LinearOp > invDiag_
(Approximate) Inverses of the diagonal operators
Teuchos::RCP< const Thyra::ProductVectorSpaceBase< double > > productRange_
Range vector space.
virtual VectorSpace range() const
Range space of this operator.