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)
81 const Teuchos::RCP<const Epetra_Operator>& content) {
82 fullContent_ = content;
83 blockedMapping_ = rcp(
new BlockedMappingStrategy(
84 vars, Teuchos::rcpFromRef(fullContent_->OperatorDomainMap()), fullContent_->Comm()));
85 SetMapStrategy(blockedMapping_);
88 BuildBlockedOperator();
91 void BlockedEpetraOperator::BuildBlockedOperator() {
92 TEUCHOS_ASSERT(blockedMapping_ != Teuchos::null);
95 const RCP<const Epetra_CrsMatrix> crsContent =
96 rcp_dynamic_cast<
const Epetra_CrsMatrix>(fullContent_);
99 if (blockedOperator_ == Teuchos::null) {
100 blockedOperator_ = blockedMapping_->buildBlockedThyraOp(crsContent, label_);
102 const RCP<Thyra::BlockedLinearOpBase<double> > blkOp =
103 rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockedOperator_,
true);
104 blockedMapping_->rebuildBlockedThyraOp(crsContent, blkOp);
108 SetOperator(blockedOperator_,
false);
111 if (reorderManager_ != Teuchos::null)
Reorder(*reorderManager_);
114 const Teuchos::RCP<const Epetra_Operator> BlockedEpetraOperator::GetBlock(
int i,
int j)
const {
115 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp =
116 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(
getThyraOp());
118 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i, j));
128 RCP<const MappingStrategy> reorderMapping =
129 rcp(
new ReorderedMappingStrategy(*reorderManager_, blockedMapping_));
130 RCP<const Thyra::BlockedLinearOpBase<double> > blockOp =
131 rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(blockedOperator_);
133 RCP<const Thyra::LinearOpBase<double> > A = buildReorderedLinearOp(*reorderManager_, blockOp);
136 SetMapStrategy(reorderMapping);
137 SetOperator(A,
false);
142 SetMapStrategy(blockedMapping_);
143 SetOperator(blockedOperator_,
false);
144 reorderManager_ = Teuchos::null;
150 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockOp =
151 rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> >(blockedOperator_);
154 int rows = Teko::blockRowCount(blockOp);
156 for (
int i = 0; i < rows; i++) {
157 for (
int j = 0; j < rows; j++) {
159 RCP<const Epetra_RowMatrix> mat = Teuchos::rcp_dynamic_cast<
const Epetra_RowMatrix>(
160 Thyra::get_Epetra_Operator(*blockOp->getBlock(i, j)));
163 EpetraExt::RowMatrixToMatrixMarketFile(formatBlockName(prefix, i, j, rows).c_str(), *mat);
169 Epetra_Vector xf(OperatorRangeMap());
170 Epetra_Vector xs(OperatorRangeMap());
171 Epetra_Vector y(OperatorDomainMap());
175 double diffNorm = 0.0, trueNorm = 0.0;
176 for (
int i = 0; i < count; i++) {
183 fullContent_->Apply(y, xf);
186 xs.Update(-1.0, xf, 1.0);
191 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