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(
73 int numVars,
const Teuchos::RCP<
const Tpetra::Operator<ST, LO, GO, NT> >& content,
74 const std::string& label)
75 : 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(
85 const std::vector<int>& vars,
86 const Teuchos::RCP<
const Tpetra::Operator<ST, LO, GO, NT> >& content,
const std::string& label)
87 : Teko::TpetraHelpers::TpetraOperatorWrapper(), label_(label) {
88 SetContent(vars, content);
91 void StridedTpetraOperator::SetContent(
92 const std::vector<int>& vars,
93 const Teuchos::RCP<
const Tpetra::Operator<ST, LO, GO, NT> >& content) {
94 fullContent_ = content;
95 stridedMapping_ = rcp(
new TpetraStridedMappingStrategy(vars, fullContent_->getDomainMap(),
96 *fullContent_->getDomainMap()->getComm()));
97 SetMapStrategy(stridedMapping_);
100 BuildBlockedOperator();
103 void StridedTpetraOperator::BuildBlockedOperator() {
104 TEUCHOS_ASSERT(stridedMapping_ != Teuchos::null);
107 const RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > crsContent =
108 rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST, LO, GO, NT> >(fullContent_);
111 if (stridedOperator_ == Teuchos::null) {
112 stridedOperator_ = stridedMapping_->buildBlockedThyraOp(crsContent, label_);
114 const RCP<Thyra::BlockedLinearOpBase<ST> > blkOp =
115 rcp_dynamic_cast<Thyra::BlockedLinearOpBase<ST> >(stridedOperator_,
true);
116 stridedMapping_->rebuildBlockedThyraOp(crsContent, blkOp);
120 SetOperator(stridedOperator_,
false);
123 if (reorderManager_ != Teuchos::null) Reorder(*reorderManager_);
126 const Teuchos::RCP<const Tpetra::Operator<ST, LO, GO, NT> > StridedTpetraOperator::GetBlock(
127 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 =
132 rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(blkOp->getBlock(i, j),
true);
133 return tOp->getConstTpetraOperator();
139 void StridedTpetraOperator::Reorder(
const BlockReorderManager& brm) {
140 reorderManager_ = rcp(
new BlockReorderManager(brm));
143 RCP<const MappingStrategy> reorderMapping =
144 rcp(
new TpetraReorderedMappingStrategy(*reorderManager_, stridedMapping_));
145 RCP<const Thyra::BlockedLinearOpBase<ST> > blockOp =
146 rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<ST> >(stridedOperator_);
148 RCP<const Thyra::LinearOpBase<ST> > A = buildReorderedLinearOp(*reorderManager_, blockOp);
151 SetMapStrategy(reorderMapping);
152 SetOperator(A,
false);
156 void StridedTpetraOperator::RemoveReording() {
157 SetMapStrategy(stridedMapping_);
158 SetOperator(stridedOperator_,
false);
159 reorderManager_ = Teuchos::null;
164 void StridedTpetraOperator::WriteBlocks(
const std::string& prefix)
const {
165 RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > blockOp =
166 rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<ST> >(stridedOperator_);
169 int rows = Teko::blockRowCount(blockOp);
171 for (
int i = 0; i < rows; i++) {
172 for (
int j = 0; j < rows; j++) {
174 RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp =
175 rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(blockOp->getBlock(i, j));
176 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > mat =
177 Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST, LO, GO, NT> >(
178 tOp->getConstTpetraOperator());
181 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<ST, LO, GO, NT> >::writeSparseFile(
182 formatBlockName(prefix, i, j, rows).c_str(), mat);
194 std::string StridedTpetraOperator::PrintNorm(
const eNormType& nrmType,
const char newline) {
195 BlockedLinearOp blockOp = toBlockedLinearOp(stridedOperator_);
198 int rowCount = Teko::blockRowCount(blockOp);
199 int colCount = Teko::blockRowCount(blockOp);
201 std::stringstream ss;
203 ss << std::scientific;
204 for (
int row = 0; row < rowCount; row++) {
205 for (
int col = 0; col < colCount; col++) {
207 RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > Aij =
208 rcp_dynamic_cast<
const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(
209 blockOp->getBlock(row, col),
true);
210 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT> > mat =
211 Teuchos::rcp_dynamic_cast<
const Tpetra::CrsMatrix<ST, LO, GO, NT> >(
212 Aij->getConstTpetraOperator(),
true);
223 case Frobenius: norm = mat->getFrobeniusNorm();
break;
225 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
226 "StridedTpetraOperator::NormType incorrectly specified. Only "
227 "Frobenius norm implemented for Tpetra matrices.");
238 bool StridedTpetraOperator::testAgainstFullOperator(
int count, ST tol)
const {
239 Tpetra::Vector<ST, LO, GO, NT> xf(getRangeMap());
240 Tpetra::Vector<ST, LO, GO, NT> xs(getRangeMap());
241 Tpetra::Vector<ST, LO, GO, NT> y(getDomainMap());
245 ST diffNorm = 0.0, trueNorm = 0.0;
246 for (
int i = 0; i < count; i++) {
253 fullContent_->apply(y, xf);
256 xs.update(-1.0, xf, 1.0);
257 diffNorm = xs.norm2();
258 trueNorm = xf.norm2();
261 result &= (diffNorm / trueNorm < tol);