47 #include "Teko_BlockedTpetraOperator.hpp"
48 #include "Teko_TpetraBlockedMappingStrategy.hpp"
49 #include "Teko_TpetraReorderedMappingStrategy.hpp"
51 #include "Teuchos_VerboseObject.hpp"
53 #include "Thyra_LinearOpBase.hpp"
54 #include "Thyra_TpetraLinearOp.hpp"
55 #include "Thyra_TpetraThyraWrappers.hpp"
56 #include "Thyra_DefaultProductMultiVector.hpp"
57 #include "Thyra_DefaultProductVectorSpace.hpp"
58 #include "Thyra_DefaultBlockedLinearOp.hpp"
60 #include "MatrixMarket_Tpetra.hpp"
65 namespace TpetraHelpers {
69 using Teuchos::rcp_dynamic_cast;
72 const std::vector<std::vector<GO> >& vars,
73 const Teuchos::RCP<
const Tpetra::Operator<ST, LO, GO, NT> >& content,
const std::string& label)
79 const std::vector<std::vector<GO> >& vars,
80 const Teuchos::RCP<
const Tpetra::Operator<ST, LO, GO, NT> >& content) {
81 fullContent_ = content;
82 blockedMapping_ = rcp(
new TpetraBlockedMappingStrategy(vars, fullContent_->getDomainMap(),
83 *fullContent_->getDomainMap()->getComm()));
84 SetMapStrategy(blockedMapping_);
87 BuildBlockedOperator();
90 void BlockedTpetraOperator::BuildBlockedOperator() {
91 TEUCHOS_ASSERT(blockedMapping_ != Teuchos::null);
94 const RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > crsContent =
95 rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST, LO, GO, NT> >(fullContent_);
98 if (blockedOperator_ == Teuchos::null) {
99 blockedOperator_ = blockedMapping_->buildBlockedThyraOp(crsContent, label_);
101 const RCP<Thyra::BlockedLinearOpBase<ST> > blkOp =
102 rcp_dynamic_cast<Thyra::BlockedLinearOpBase<ST> >(blockedOperator_,
true);
103 blockedMapping_->rebuildBlockedThyraOp(crsContent, blkOp);
107 SetOperator(blockedOperator_,
false);
110 if (reorderManager_ != Teuchos::null)
Reorder(*reorderManager_);
113 const Teuchos::RCP<const Tpetra::Operator<ST, LO, GO, NT> > BlockedTpetraOperator::GetBlock(
114 int i,
int j)
const {
115 const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp =
116 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<ST> >(
getThyraOp());
118 RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp =
119 rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(blkOp->getBlock(i, j),
true);
120 return tOp->getConstTpetraOperator();
130 RCP<const MappingStrategy> reorderMapping =
131 rcp(
new TpetraReorderedMappingStrategy(*reorderManager_, blockedMapping_));
132 RCP<const Thyra::BlockedLinearOpBase<ST> > blockOp =
133 rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<ST> >(blockedOperator_);
135 RCP<const Thyra::LinearOpBase<ST> > A = buildReorderedLinearOp(*reorderManager_, blockOp);
138 SetMapStrategy(reorderMapping);
139 SetOperator(A,
false);
144 SetMapStrategy(blockedMapping_);
145 SetOperator(blockedOperator_,
false);
146 reorderManager_ = Teuchos::null;
152 RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > blockOp =
153 rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<ST> >(blockedOperator_);
156 int rows = Teko::blockRowCount(blockOp);
158 for (
int i = 0; i < rows; i++) {
159 for (
int j = 0; j < rows; j++) {
161 RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp =
162 rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(blockOp->getBlock(i, j));
163 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > mat =
164 Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST, LO, GO, NT> >(
165 tOp->getConstTpetraOperator());
168 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<ST, LO, GO, NT> >::writeSparseFile(
169 formatBlockName(prefix, i, j, rows).c_str(), mat);
175 Tpetra::Vector<ST, LO, GO, NT> xf(getRangeMap());
176 Tpetra::Vector<ST, LO, GO, NT> xs(getRangeMap());
177 Tpetra::Vector<ST, LO, GO, NT> y(getDomainMap());
181 ST diffNorm = 0.0, trueNorm = 0.0;
182 for (
int i = 0; i < count; i++) {
189 fullContent_->apply(y, xf);
192 xs.update(-1.0, xf, 1.0);
193 diffNorm = xs.norm2();
194 trueNorm = xf.norm2();
197 result &= (diffNorm / trueNorm < tol);
void Reorder(const BlockReorderManager &brm)
Class that describes how a flat blocked operator should be reordered.
void RemoveReording()
Remove any reordering on this object.
const RCP< const Thyra::LinearOpBase< ST > > getThyraOp() const
Return the thyra operator associated with this wrapper.
virtual void WriteBlocks(const std::string &prefix) const
BlockedTpetraOperator(const std::vector< std::vector< GO > > &vars, const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &content, const std::string &label="<ANYM>")
bool testAgainstFullOperator(int count, ST tol) const
Helps perform sanity checks.
virtual void SetContent(const std::vector< std::vector< GO > > &vars, const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &content)
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...