47 #include "Teko_BlockedEpetraOperator.hpp" 48 #include "Teko_BlockedMappingStrategy.hpp" 49 #include "Teko_ReorderedMappingStrategy.hpp" 51 #include "Teuchos_VerboseObject.hpp" 53 #include "Thyra_LinearOpBase.hpp" 54 #include "Thyra_EpetraLinearOp.hpp" 55 #include "Thyra_EpetraThyraWrappers.hpp" 56 #include "Thyra_DefaultProductMultiVector.hpp" 57 #include "Thyra_DefaultProductVectorSpace.hpp" 58 #include "Thyra_DefaultBlockedLinearOp.hpp" 59 #include "Thyra_get_Epetra_Operator.hpp" 61 #include "EpetraExt_MultiVectorOut.h" 62 #include "EpetraExt_RowMatrixOut.h" 71 using Teuchos::rcp_dynamic_cast;
74 const Teuchos::RCP<const Epetra_Operator> & content,
75 const std::string & label)
82 const Teuchos::RCP<const Epetra_Operator> & content)
84 fullContent_ = content;
85 blockedMapping_ = rcp(
new BlockedMappingStrategy(vars,Teuchos::rcpFromRef(fullContent_->OperatorDomainMap()),
86 fullContent_->Comm()));
87 SetMapStrategy(blockedMapping_);
90 BuildBlockedOperator();
93 void BlockedEpetraOperator::BuildBlockedOperator()
95 TEUCHOS_ASSERT(blockedMapping_!=Teuchos::null);
98 const RCP<const Epetra_CrsMatrix> crsContent
99 = rcp_dynamic_cast<
const Epetra_CrsMatrix>(fullContent_);
102 if(blockedOperator_==Teuchos::null) {
103 blockedOperator_ = blockedMapping_->buildBlockedThyraOp(crsContent,label_);
106 const RCP<Thyra::BlockedLinearOpBase<double> > blkOp
107 = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockedOperator_,
true);
108 blockedMapping_->rebuildBlockedThyraOp(crsContent,blkOp);
112 SetOperator(blockedOperator_,
false);
115 if(reorderManager_!=Teuchos::null)
119 const Teuchos::RCP<const Epetra_Operator> BlockedEpetraOperator::GetBlock(
int i,
int j)
const 121 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
122 = Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(
getThyraOp());
124 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i,j));
135 RCP<const MappingStrategy> reorderMapping = rcp(
new ReorderedMappingStrategy(*reorderManager_,blockedMapping_));
136 RCP<const Thyra::BlockedLinearOpBase<double> > blockOp
137 = rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(blockedOperator_);
139 RCP<const Thyra::LinearOpBase<double> > A = buildReorderedLinearOp(*reorderManager_,blockOp);
142 SetMapStrategy(reorderMapping);
143 SetOperator(A,
false);
149 SetMapStrategy(blockedMapping_);
150 SetOperator(blockedOperator_,
false);
151 reorderManager_ = Teuchos::null;
158 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockOp
159 = rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> >(blockedOperator_);
164 for(
int i=0;i<rows;i++) {
165 for(
int j=0;j<rows;j++) {
167 std::stringstream ss;
168 ss << prefix <<
"_" << i << j <<
".mm";
171 RCP<const Epetra_RowMatrix> mat
172 = Teuchos::rcp_dynamic_cast<
const Epetra_RowMatrix>(Thyra::get_Epetra_Operator(*blockOp->getBlock(i,j)));
175 EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*mat);
180 #ifndef Teko_DEBUG_OFF 183 Epetra_Vector xf(OperatorRangeMap());
184 Epetra_Vector xs(OperatorRangeMap());
185 Epetra_Vector y(OperatorDomainMap());
189 double 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);
205 result &= (diffNorm/trueNorm < tol);
const RCP< const Thyra::LinearOpBase< double > > getThyraOp() const
Return the thyra operator associated with this wrapper.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
virtual void WriteBlocks(const std::string &prefix) const
bool testAgainstFullOperator(int count, double tol) const
Helps perform sanity checks.
Class that describes how a flat blocked operator should be reordered.
void Reorder(const BlockReorderManager &brm)
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...
virtual void SetContent(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< const Epetra_Operator > &content)
void RemoveReording()
Remove any reordering on this object.
BlockedEpetraOperator(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< const Epetra_Operator > &content, const std::string &label="<ANYM>")