Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_ThyraUtils_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Xpetra: A linear algebra interface package
4//
5// Copyright 2012 NTESS and the Xpetra contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef HAVE_XPETRA_THYRAUTILS_DEF_HPP
11#define HAVE_XPETRA_THYRAUTILS_DEF_HPP
12
13#ifdef HAVE_XPETRA_THYRA
14
15#include "Xpetra_BlockedCrsMatrix.hpp"
16
17#include "Xpetra_ThyraUtils_decl.hpp"
18
19namespace Xpetra {
20
21template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
22Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>
23Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
24 toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm, std::vector<size_t>& stridingInfo, LocalOrdinal stridedBlockId, GlobalOrdinal offset) {
25 Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map = toXpetra(vectorSpace, comm);
26
27 if (stridedBlockId == -1) {
28 TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo.size() != 0);
29 } else {
30 TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo[stridedBlockId] != 0);
31 }
32
33 Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>> ret = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(map, stridingInfo, stridedBlockId, offset);
34 return ret;
35}
36
37template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
38Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
39Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
40 toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
41 using Teuchos::as;
42 using Teuchos::RCP;
43 using Teuchos::rcp_dynamic_cast;
44 using ThyVecSpaceBase = Thyra::VectorSpaceBase<Scalar>;
45 using ThyProdVecSpaceBase = Thyra::ProductVectorSpaceBase<Scalar>;
46 using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
47
48 RCP<Map> resultMap = Teuchos::null;
49 RCP<const ThyProdVecSpaceBase> prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase>(vectorSpace);
50 if (prodVectorSpace != Teuchos::null) {
51 // SPECIAL CASE: product Vector space
52 // collect all submaps to store them in a hierarchical BlockedMap object
53 TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks() == 0, std::logic_error, "Found a product vector space with zero blocks.");
54 std::vector<RCP<const Map>> mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
55 std::vector<RCP<const Map>> mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
56 for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
57 RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
58 // can be of type Map or BlockedMap (containing Thyra GIDs)
59 mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
60 }
61
62 // get offsets for submap GIDs
63 // we need that for the full map (Xpetra GIDs)
64 std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(), 0);
65 for (int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
66 gidOffsets[i] = mapsThyra[i - 1]->getMaxAllGlobalIndex() + gidOffsets[i - 1] + 1;
67 }
68
69 for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
70 RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
71 // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
72 mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
73 }
74
75 resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapsXpetra, mapsThyra));
76 } else {
77 // STANDARD CASE: no product map
78 // check whether we have a Tpetra based Thyra operator
79 Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vectorSpace);
80 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(tpetra_vsc) == true, Xpetra::Exceptions::RuntimeError, "toXpetra failed to convert provided vector space to Thyra::TpetraVectorSpace. This is the general implementation for Tpetra only.");
81
82 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> TpMap;
83 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpVector;
84 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
85 typedef Thyra::VectorBase<Scalar> ThyVecBase;
86 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
87 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
88 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
89 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
90 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
91 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
92 resultMap = Xpetra::toXpetraNonConst(rgTpetraMap);
93 } // end standard case (no product map)
94 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
95 return resultMap;
96}
97
98template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
99Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
100Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
101 toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
102 using Teuchos::as;
103 using Teuchos::RCP;
104 using Teuchos::rcp_dynamic_cast;
105
106 using ThyMultVecBase = Thyra::MultiVectorBase<Scalar>;
107 using ThyProdMultVecBase = Thyra::ProductMultiVectorBase<Scalar>;
108 using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
109
110 // return value
111 RCP<MultiVector> xpMultVec = Teuchos::null;
112
113 // check whether v is a product multi vector
114 Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase>(v);
115 if (thyProdVec != Teuchos::null) {
116 // SPECIAL CASE: create a nested BlockedMultiVector
117 // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
118 RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
119 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
120
121 // create new Xpetra::BlockedMultiVector
122 xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
123
124 RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec, true);
125
126 // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
127 for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b) {
128 RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
129 // xpBlockMV can be of type MultiVector or BlockedMultiVector
130 RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); // recursive call
131 xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
132 }
133 } else {
134 // STANDARD CASE: no product vector
135 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> ConverterT;
136 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpMultVec;
137 typedef Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpTpMultVec;
138 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
139 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
140
141 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
142 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
143 TEUCHOS_TEST_FOR_EXCEPTION(thyraTpetraMultiVector == Teuchos::null, Xpetra::Exceptions::RuntimeError, "toXpetra failed to convert provided multi vector to Thyra::TpetraMultiVector. This is the general implementation for Tpetra only.");
144 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
145 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
146 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
147 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
148 xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
149 } // end standard case
150 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
151 return xpMultVec;
152}
153
154template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
155Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
156Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
157 toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
158 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> cv =
159 Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar>>(v);
160 Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> r =
161 toXpetra(cv, comm);
162 return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(r);
163}
164
165template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
166bool Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
167 isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
168 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(op));
169
170 // check whether we have a Tpetra based Thyra operator
171 bool bIsTpetra = false;
172 Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(op);
173 bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
174
175 // for debugging purposes: find out why dynamic cast failed
176 if (!bIsTpetra &&
177 Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar>>(op) == Teuchos::null) {
178 // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
179 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
180 if (Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
181 std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
182 std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
183 std::cout << " properly set!" << std::endl;
184 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
185 }
186 }
187
188#if 0
189 // Check whether it is a blocked operator.
190 // If yes, grab the (0,0) block and check the underlying linear algebra
191 // Note: we expect that the (0,0) block exists!
192 if(bIsTpetra == false) {
193 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
194 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
195 if(ThyBlockedOp != Teuchos::null) {
196 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
197 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
198 ThyBlockedOp->getBlock(0,0);
199 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
200 bIsTpetra = isTpetra(b00);
201 }
202 }
203#endif
204
205 return bIsTpetra;
206}
207
208template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
209bool Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
210 isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
211 // Check whether it is a blocked operator.
212 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar>> ThyBlockedOp =
213 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar>>(op);
214 if (ThyBlockedOp != Teuchos::null) {
215 return true;
216 }
217 return false;
218}
219
220template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
221Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
222Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
223 toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
224 if (isTpetra(op)) {
225 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
226 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
227 // we should also add support for the const versions!
228 // getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
229 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
230 Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
231 Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
232 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat);
233 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
234
235 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
236 Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraNcnstCrsMat));
237 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
238
239 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
240 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
241 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
242 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
243 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
244 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
245 return xpMat;
246 }
247
248 return Teuchos::null;
249}
250
251template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
252Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
253Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
254 toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
255 using Teuchos::RCP;
256 using Teuchos::rcp;
257 using Teuchos::rcp_const_cast;
258 using Teuchos::rcp_dynamic_cast;
259
260 if (isTpetra(op)) {
261 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
262 typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraOperator_t;
263 RCP<const TpetraOperator_t> TpetraOp = TOE::getConstTpetraOperator(op);
264 TEUCHOS_TEST_FOR_EXCEPT(is_null(TpetraOp));
265 RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat =
266 rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
267 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat =
268 rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
269
270 RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
271 rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(
272 rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat)));
273 TEUCHOS_TEST_FOR_EXCEPT(is_null(xTpetraCrsMat));
274 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
275 rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
276 RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
277 rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
278 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
279 rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
280 return xpMat;
281 }
282
283 return Teuchos::null;
284}
285
286template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
287Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
288Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
289 toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
290 if (isTpetra(op)) {
291 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
292 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
293
294 auto nonConstTpetraOp = Teuchos::rcp_const_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp);
295
296 Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
297 Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(nonConstTpetraOp));
298 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
299
300 Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
301 Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
302 return xpOp;
303 }
304
305 return Teuchos::null;
306}
307
308template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
309Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
310Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
311 toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
312 if (isTpetra(op)) {
313 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
314 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
315
316 Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
317 Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
318 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
319
320 Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
321 Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
322 return xpOp;
323 }
324
325 return Teuchos::null;
326}
327
328template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
329Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
330Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
331 toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op) {
332 using Teuchos::rcp_const_cast;
333 using Teuchos::rcp_dynamic_cast;
334
335 RCP<const Thyra::VectorBase<Scalar>> diag = op->getDiag();
336
337 RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpDiag;
338 using thyTpV = Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
339 using tV = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
340 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
341 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
342 if (!tDiag.is_null())
343 xpDiag = Xpetra::toXpetra(tDiag);
344 }
345 TEUCHOS_ASSERT(!xpDiag.is_null());
346 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
347 return M;
348}
349
350template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
351Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
352Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
353 toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op) {
354 return toXpetra(Teuchos::rcp_const_cast<Thyra::DiagonalLinearOpBase<Scalar>>(op));
355}
356
357template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
358Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
359Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
360 toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map) {
361 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> thyraMap = Teuchos::null;
362
363 // check whether map is of type BlockedMap
364 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
365 if (bmap.is_null() == false) {
366 Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>> vecSpaces(bmap->getNumMaps());
367 for (size_t i = 0; i < bmap->getNumMaps(); i++) {
368 // we need Thyra GIDs for all the submaps
369 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> vs =
370 Xpetra::ThyraUtils<Scalar, LO, GO, Node>::toThyra(bmap->getMap(i, true));
371 vecSpaces[i] = vs;
372 }
373
374 thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
375 return thyraMap;
376 }
377
378 // standard case
379 if (map->lib() == Xpetra::UseTpetra) {
380 Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>> tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(map);
381 if (tpetraMap == Teuchos::null)
382 throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
383 RCP<Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetraMap->getTpetra_Map());
384 thyraMap = thyraTpetraMap;
385 }
386
387 return thyraMap;
388}
389
390template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
391Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
392Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
393 toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
394 // create Thyra MultiVector
395 if (vec->getMap()->lib() == Xpetra::UseTpetra) {
396 auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
397 RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpMV = Teuchos::rcp_dynamic_cast<const TpetraMultiVector>(vec)->getTpetra_MultiVector();
398 auto thyDomMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(vec->getNumVectors(), vec->getMap()->getComm()));
399 auto thyMV = rcp(new Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
400 thyMV->initialize(thyTpMap, thyDomMap, tpMV);
401 return thyMV;
402 }
403
404 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "MultiVector cannot be converted to Thyra.");
405}
406
407template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
408Teuchos::RCP<const Thyra::VectorBase<Scalar>>
409Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
410 toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
411 // create Thyra Vector
412 if (vec->getMap()->lib() == Xpetra::UseTpetra) {
413 auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
414 RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpVec = Teuchos::rcp_dynamic_cast<const TpetraVector>(vec)->getTpetra_Vector();
415 auto thyVec = rcp(new Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
416 thyVec->initialize(thyTpMap, tpVec);
417 return thyVec;
418 }
419
420 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Vector cannot be converted to Thyra.");
421}
422
423template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
424void Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
425 updateThyra(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> source, Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node>> mapExtractor, const Teuchos::RCP<Thyra::MultiVectorBase<Scalar>>& target) {
426 using Teuchos::as;
427 using Teuchos::RCP;
428 using Teuchos::rcp_dynamic_cast;
429 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
430 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
431 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
432 // typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
433 // typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
434 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
435
436 // copy data from tY_inout to Y_inout
437 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
438 if (prodTarget != Teuchos::null) {
439 RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
440 if (bSourceVec.is_null() == true) {
441 // SPECIAL CASE: target vector is product vector:
442 // update Thyra product multi vector with data from a merged Xpetra multi vector
443
444 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
445 TEUCHOS_TEST_FOR_EXCEPTION(prodTarget->productSpace()->numBlocks() != as<int>(mapExtractor->NumMaps()), std::logic_error, "Inconsistent numbers of sub maps in Thyra::ProductVectorSpace and Xpetra::MapExtractor.");
446
447 for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
448 // access Xpetra data
449 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
450
451 // access Thyra data
452 Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
453 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
454 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
455 const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
456 const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim());
457 RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
458 Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
459
460 // loop over all vectors in multivector
461 for (size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
462 Teuchos::ArrayRCP<const Scalar> xpData = xpSubBlock->getData(j); // access const data from Xpetra object
463
464 // loop over all local rows
465 for (LocalOrdinal i = 0; i < localSubDim; ++i) {
466 (*thyData)(i, j) = xpData[i];
467 }
468 }
469 }
470 } else {
471 // source vector is a blocked multivector
472 // TODO test me
473 TEUCHOS_TEST_FOR_EXCEPTION(prodTarget->productSpace()->numBlocks() != as<int>(bSourceVec->getBlockedMap()->getNumMaps()), std::logic_error, "Inconsistent numbers of sub maps in Thyra::ProductVectorSpace and Xpetra::BlockedMultiVector.");
474
475 for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
476 // access Thyra data
477 RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
478
479 Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
480
481 // access Thyra data
482 Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
483 Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
484 }
485 }
486 } else {
487 // STANDARD case:
488 // update Thyra::MultiVector with data from an Xpetra::MultiVector
489
490 // access Thyra data
491 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
492 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
493 const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
494 const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim());
495 RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
496 Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
497
498 // loop over all vectors in multivector
499 for (size_t j = 0; j < source->getNumVectors(); ++j) {
500 Teuchos::ArrayRCP<const Scalar> xpData = source->getData(j); // access const data from Xpetra object
501 // loop over all local rows
502 for (LocalOrdinal i = 0; i < localSubDim; ++i) {
503 (*thyData)(i, j) = xpData[i];
504 }
505 }
506 }
507}
508
509template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
510Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
511Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
512 toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
513 // create a Thyra operator from Xpetra::CrsMatrix
514 Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
515
516 // bool bIsTpetra = false;
517
518 Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
519 if (tpetraMat != Teuchos::null) {
520 // bIsTpetra = true;
521 Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
522 Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
523 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
524
525 Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
526 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
527
528 thyraOp = Thyra::createConstLinearOp(tpOperator);
529 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
530 } else {
531 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to Tpetra::CrsMatrix failed. Assume matrix should be Epetra then. No Epetra available");
532 }
533 return thyraOp;
534}
535
536template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
537Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
538Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
539 toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
540 // create a Thyra operator from Xpetra::CrsMatrix
541 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
542
543 // bool bIsTpetra = false;
544
545 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
546 if (tpetraMat != Teuchos::null) {
547 // bIsTpetra = true;
548 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
549 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
550 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
551
552 Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
553 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
554
555 thyraOp = Thyra::createLinearOp(tpOperator);
556 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
557 } else {
558 // cast to TpetraCrsMatrix failed
559 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Guess, matrix should be Epetra then, but no Epetra available.");
560 }
561 return thyraOp;
562}
563
564template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
565Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
566Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
567 toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
568 int nRows = mat->Rows();
569 int nCols = mat->Cols();
570
571 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Ablock = mat->getInnermostCrsMatrix();
572 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock);
573 TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
574
575 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
576 if (tpetraMat != Teuchos::null) {
577 // create new Thyra blocked operator
578 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar>> blockMat =
579 Thyra::defaultBlockedLinearOp<Scalar>();
580
581 blockMat->beginBlockFill(nRows, nCols);
582
583 for (int r = 0; r < nRows; ++r) {
584 for (int c = 0; c < nCols; ++c) {
585 Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r, c);
586
587 if (xpmat == Teuchos::null) continue; // shortcut for empty blocks
588
589 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thBlock = Teuchos::null;
590
591 // check whether the subblock is again a blocked operator
592 Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpblock =
593 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpmat);
594 if (xpblock != Teuchos::null) {
595 if (xpblock->Rows() == 1 && xpblock->Cols() == 1) {
596 // If it is a single block operator, unwrap it
597 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
598 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
599 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
600 } else {
601 // recursive call for general blocked operators
602 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpblock);
603 }
604 } else {
605 // check whether it is a CRSMatrix object
606 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
607 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
608 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
609 }
610
611 blockMat->setBlock(r, c, thBlock);
612 }
613 }
614
615 blockMat->endBlockFill();
616
617 return blockMat;
618 } else {
619 // tpetraMat == Teuchos::null
620 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Guess, matrix should be Epetra then, but no Epetra available.");
621 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
622 }
623}
624
625} // namespace Xpetra
626
627#endif
628
629#endif
Xpetra_ThyraUtils_decl.hpp
Teuchos::RCP::is_null
bool is_null() const
Teuchos::RCP::ptr
Ptr< T > ptr() const
Xpetra::Exceptions::BadCast
Exception indicating invalid cast attempted.
Definition Xpetra_Exceptions.hpp:42
Xpetra::MapUtils::transformThyra2XpetraGIDs
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > transformThyra2XpetraGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlInput, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlReferenceInput)
replace set of global ids by new global ids
Definition Xpetra_MapUtils.hpp:191
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
Definition Xpetra_MultiVectorFactory_def.hpp:24
Xpetra::StridedMapFactory::Build
static RCP< Xpetra::StridedMap< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
Definition Xpetra_StridedMapFactory_def.hpp:24
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
Teuchos::as
TypeTo as(const TypeFrom &t)
TEUCHOS_UNREACHABLE_RETURN
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Tpetra::createLocalMapWithNode
Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > createLocalMapWithNode(const size_t numElements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm)
Xpetra::toXpetraNonConst
const RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > toXpetraNonConst(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &map)
Definition Xpetra_TpetraMap_decl.hpp:224
Xpetra::toXpetra
RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toXpetra(RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > graph)
Definition Xpetra_TpetraCrsGraph_decl.hpp:365
Xpetra::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)