47 #include "Teko_StridedTpetraOperator.hpp"
48 #include "Teko_TpetraStridedMappingStrategy.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 "Tpetra_Vector.hpp"
61 #include "MatrixMarket_Tpetra.hpp"
66 namespace TpetraHelpers {
70 using Teuchos::rcp_dynamic_cast;
72 StridedTpetraOperator::StridedTpetraOperator(
int numVars,
const Teuchos::RCP<
const Tpetra::Operator<ST,LO,GO,NT> > & content,
73 const std::string & label)
74 : Teko::TpetraHelpers::TpetraOperatorWrapper(), label_(label)
76 std::vector<int> vars;
79 for(
int i=0;i<numVars;i++) vars.push_back(1);
81 SetContent(vars,content);
84 StridedTpetraOperator::StridedTpetraOperator(
const std::vector<int> & vars,
const Teuchos::RCP<
const Tpetra::Operator<ST,LO,GO,NT> > & content,
85 const std::string & label)
86 : Teko::TpetraHelpers::TpetraOperatorWrapper(), label_(label)
88 SetContent(vars,content);
91 void StridedTpetraOperator::SetContent(
const std::vector<int> & vars,
const Teuchos::RCP<
const Tpetra::Operator<ST,LO,GO,NT> > & content)
93 fullContent_ = content;
94 stridedMapping_ = rcp(
new TpetraStridedMappingStrategy(vars,fullContent_->getDomainMap(),
95 *fullContent_->getDomainMap()->getComm()));
96 SetMapStrategy(stridedMapping_);
99 BuildBlockedOperator();
102 void StridedTpetraOperator::BuildBlockedOperator()
104 TEUCHOS_ASSERT(stridedMapping_!=Teuchos::null);
107 const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsContent = rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(fullContent_);
110 if(stridedOperator_==Teuchos::null) {
111 stridedOperator_ = stridedMapping_->buildBlockedThyraOp(crsContent,label_);
114 const RCP<Thyra::BlockedLinearOpBase<ST> > blkOp = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<ST> >(stridedOperator_,
true);
115 stridedMapping_->rebuildBlockedThyraOp(crsContent,blkOp);
119 SetOperator(stridedOperator_,
false);
122 if(reorderManager_!=Teuchos::null)
123 Reorder(*reorderManager_);
126 const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > StridedTpetraOperator::GetBlock(
int i,
int j)
const
128 const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
129 = Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());
131 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blkOp->getBlock(i,j),
true);
132 return tOp->getConstTpetraOperator();
138 void StridedTpetraOperator::Reorder(
const BlockReorderManager & brm)
140 reorderManager_ = rcp(
new BlockReorderManager(brm));
143 RCP<const MappingStrategy> reorderMapping = rcp(
new TpetraReorderedMappingStrategy(*reorderManager_,stridedMapping_));
144 RCP<const Thyra::BlockedLinearOpBase<ST> > blockOp
145 = rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<ST> >(stridedOperator_);
147 RCP<const Thyra::LinearOpBase<ST> > A = buildReorderedLinearOp(*reorderManager_,blockOp);
150 SetMapStrategy(reorderMapping);
151 SetOperator(A,
false);
155 void StridedTpetraOperator::RemoveReording()
157 SetMapStrategy(stridedMapping_);
158 SetOperator(stridedOperator_,
false);
159 reorderManager_ = Teuchos::null;
164 void StridedTpetraOperator::WriteBlocks(
const std::string & prefix)
const
166 RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > blockOp
167 = rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<ST> >(stridedOperator_);
170 int rows = Teko::blockRowCount(blockOp);
172 for(
int i=0;i<rows;i++) {
173 for(
int j=0;j<rows;j++) {
175 std::stringstream ss;
176 ss << prefix <<
"_" << i << j <<
".mm";
179 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blockOp->getBlock(i,j));
180 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > mat
181 = Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator());
184 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<ST,LO,GO,NT> >::writeSparseFile(ss.str().c_str(),mat);
196 std::string StridedTpetraOperator::PrintNorm(
const eNormType & nrmType,
const char newline)
198 BlockedLinearOp blockOp = toBlockedLinearOp(stridedOperator_);
201 int rowCount = Teko::blockRowCount(blockOp);
202 int colCount = Teko::blockRowCount(blockOp);
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 Thyra::TpetraLinearOp<ST,LO,GO,NT> > Aij = rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blockOp->getBlock(row,col),
true);
211 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > mat = Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST,LO,GO,NT> >(Aij->getConstTpetraOperator(),
true);
223 norm = mat->getFrobeniusNorm();
226 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::runtime_error,
227 "StridedTpetraOperator::NormType incorrectly specified. Only Frobenius norm implemented for Tpetra matrices.");
238 #ifndef Teko_DEBUG_OFF
239 bool StridedTpetraOperator::testAgainstFullOperator(
int count,ST tol)
const
241 Tpetra::Vector<ST,LO,GO,NT> xf(getRangeMap());
242 Tpetra::Vector<ST,LO,GO,NT> xs(getRangeMap());
243 Tpetra::Vector<ST,LO,GO,NT> y(getDomainMap());
247 ST 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);
259 diffNorm = xs.norm2();
260 trueNorm = xf.norm2();
263 result &= (diffNorm/trueNorm < tol);