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_);
162 int rows = Teko::blockRowCount(blockOp);
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);
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>")
const RCP< const Thyra::LinearOpBase< double > > getThyraOp() const
Return the thyra operator associated with this wrapper.
virtual void WriteBlocks(const std::string &prefix) const