47 #include "Epetra/Teko_StridedEpetraOperator.hpp" 48 #include "Epetra/Teko_StridedMappingStrategy.hpp" 49 #include "Epetra/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 "Epetra_Vector.h" 62 #include "EpetraExt_MultiVectorOut.h" 63 #include "EpetraExt_RowMatrixOut.h" 72 using Teuchos::rcp_dynamic_cast;
74 StridedEpetraOperator::StridedEpetraOperator(
int numVars,
const Teuchos::RCP<const Epetra_Operator> & content,
75 const std::string & label)
76 :
Teko::Epetra::EpetraOperatorWrapper(), label_(label)
78 std::vector<int> vars;
81 for(
int i=0;i<numVars;i++) vars.push_back(1);
83 SetContent(vars,content);
86 StridedEpetraOperator::StridedEpetraOperator(
const std::vector<int> & vars,
const Teuchos::RCP<const Epetra_Operator> & content,
87 const std::string & label)
88 :
Teko::Epetra::EpetraOperatorWrapper(), label_(label)
90 SetContent(vars,content);
93 void StridedEpetraOperator::SetContent(
const std::vector<int> & vars,
const Teuchos::RCP<const Epetra_Operator> & content)
95 fullContent_ = content;
96 stridedMapping_ = rcp(
new StridedMappingStrategy(vars,Teuchos::rcpFromRef(fullContent_->OperatorDomainMap()),
97 fullContent_->Comm()));
98 SetMapStrategy(stridedMapping_);
101 BuildBlockedOperator();
104 void StridedEpetraOperator::BuildBlockedOperator()
106 TEUCHOS_ASSERT(stridedMapping_!=Teuchos::null);
109 const RCP<const Epetra_CrsMatrix> crsContent = rcp_dynamic_cast<
const Epetra_CrsMatrix>(fullContent_);
112 if(stridedOperator_==Teuchos::null) {
113 stridedOperator_ = stridedMapping_->buildBlockedThyraOp(crsContent,label_);
116 const RCP<Thyra::BlockedLinearOpBase<double> > blkOp = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(stridedOperator_,
true);
117 stridedMapping_->rebuildBlockedThyraOp(crsContent,blkOp);
121 SetOperator(stridedOperator_,
false);
124 if(reorderManager_!=Teuchos::null)
125 Reorder(*reorderManager_);
128 const Teuchos::RCP<const Epetra_Operator> StridedEpetraOperator::GetBlock(
int i,
int j)
const 130 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
131 = Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
133 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i,j));
139 void StridedEpetraOperator::Reorder(
const BlockReorderManager & brm)
141 reorderManager_ = rcp(
new BlockReorderManager(brm));
144 RCP<const MappingStrategy> reorderMapping = rcp(
new ReorderedMappingStrategy(*reorderManager_,stridedMapping_));
145 RCP<const Thyra::BlockedLinearOpBase<double> > blockOp
146 = rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(stridedOperator_);
148 RCP<const Thyra::LinearOpBase<double> > A = buildReorderedLinearOp(*reorderManager_,blockOp);
151 SetMapStrategy(reorderMapping);
152 SetOperator(A,
false);
156 void StridedEpetraOperator::RemoveReording()
158 SetMapStrategy(stridedMapping_);
159 SetOperator(stridedOperator_,
false);
160 reorderManager_ = Teuchos::null;
165 void StridedEpetraOperator::WriteBlocks(
const std::string & prefix)
const 167 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockOp
168 = rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> >(stridedOperator_);
173 for(
int i=0;i<rows;i++) {
174 for(
int j=0;j<rows;j++) {
176 std::stringstream ss;
177 ss << prefix <<
"_" << i << j <<
".mm";
180 RCP<const Epetra_RowMatrix> mat
181 = Teuchos::rcp_dynamic_cast<
const Epetra_RowMatrix>(Thyra::get_Epetra_Operator(*blockOp->getBlock(i,j)));
184 EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*mat);
196 std::string StridedEpetraOperator::PrintNorm(
const eNormType & nrmType,
const char newline)
204 std::stringstream ss;
206 ss << std::scientific;
207 for(
int row=0;row<rowCount;row++) {
208 for(
int col=0;col<colCount;col++) {
210 RCP<const Epetra_CrsMatrix> mat = Teuchos::rcp_dynamic_cast<
const Epetra_CrsMatrix>(
211 Thyra::get_Epetra_Operator(*blockOp->getBlock(row,col)));
217 norm = mat->NormInf();
220 norm = mat->NormOne();
223 norm = mat->NormFrobenius();
226 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error,
227 "StridedEpetraOperator::eNormType incorrectly specified");
238 #ifndef Teko_DEBUG_OFF 239 bool StridedEpetraOperator::testAgainstFullOperator(
int count,
double tol)
const 241 Epetra_Vector xf(OperatorRangeMap());
242 Epetra_Vector xs(OperatorRangeMap());
243 Epetra_Vector y(OperatorDomainMap());
247 double diffNorm=0.0,trueNorm=0.0;
248 for(
int i=0;i<count;i++) {
255 fullContent_->Apply(y,xf);
258 xs.Update(-1.0,xf,1.0);
263 result &= (diffNorm/trueNorm < tol);
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
BlockedLinearOp toBlockedLinearOp(LinearOp &clo)
Converts a LinearOp to a BlockedLinearOp.