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 Teuchos::RCP<
const Tpetra::Operator<ST,LO,GO,NT> > & content,
73 const std::string & label)
80 const Teuchos::RCP<
const Tpetra::Operator<ST,LO,GO,NT> > & content)
82 fullContent_ = content;
83 blockedMapping_ = rcp(
new TpetraBlockedMappingStrategy(vars,fullContent_->getDomainMap(),
84 *fullContent_->getDomainMap()->getComm()));
85 SetMapStrategy(blockedMapping_);
88 BuildBlockedOperator();
91 void BlockedTpetraOperator::BuildBlockedOperator()
93 TEUCHOS_ASSERT(blockedMapping_!=Teuchos::null);
96 const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsContent
97 = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(fullContent_);
100 if(blockedOperator_==Teuchos::null) {
101 blockedOperator_ = blockedMapping_->buildBlockedThyraOp(crsContent,label_);
104 const RCP<Thyra::BlockedLinearOpBase<ST> > blkOp
105 = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<ST> >(blockedOperator_,
true);
106 blockedMapping_->rebuildBlockedThyraOp(crsContent,blkOp);
110 SetOperator(blockedOperator_,
false);
113 if(reorderManager_!=Teuchos::null)
117 const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > BlockedTpetraOperator::GetBlock(
int i,
int j)
const
119 const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
120 = Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<ST> >(
getThyraOp());
122 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blkOp->getBlock(i,j),
true);
123 return tOp->getConstTpetraOperator();
134 RCP<const MappingStrategy> reorderMapping = rcp(
new TpetraReorderedMappingStrategy(*reorderManager_,blockedMapping_));
135 RCP<const Thyra::BlockedLinearOpBase<ST> > blockOp
136 = rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<ST> >(blockedOperator_);
138 RCP<const Thyra::LinearOpBase<ST> > A = buildReorderedLinearOp(*reorderManager_,blockOp);
141 SetMapStrategy(reorderMapping);
142 SetOperator(A,
false);
148 SetMapStrategy(blockedMapping_);
149 SetOperator(blockedOperator_,
false);
150 reorderManager_ = Teuchos::null;
157 RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > blockOp
158 = rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<ST> >(blockedOperator_);
161 int rows = Teko::blockRowCount(blockOp);
163 for(
int i=0;i<rows;i++) {
164 for(
int j=0;j<rows;j++) {
166 std::stringstream ss;
167 ss << prefix <<
"_" << i << j <<
".mm";
170 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blockOp->getBlock(i,j));
171 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > mat
172 = Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator());
175 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<ST,LO,GO,NT> >::writeSparseFile(ss.str().c_str(),mat);
180 #ifndef Teko_DEBUG_OFF
183 Tpetra::Vector<ST,LO,GO,NT> xf(getRangeMap());
184 Tpetra::Vector<ST,LO,GO,NT> xs(getRangeMap());
185 Tpetra::Vector<ST,LO,GO,NT> y(getDomainMap());
189 ST diffNorm=0.0,trueNorm=0.0;
190 for(
int i=0;i<count;i++) {
197 fullContent_->apply(y,xf);
200 xs.update(-1.0,xf,1.0);
201 diffNorm = xs.norm2();
202 trueNorm = xf.norm2();
205 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...