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#ifdef HAVE_XPETRA_TPETRA
78 // STANDARD CASE: no product map
79 // check whether we have a Tpetra based Thyra operator
80 Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vectorSpace);
81 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.");
82
83 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> TpMap;
84 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpVector;
85 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
86 typedef Thyra::VectorBase<Scalar> ThyVecBase;
87 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
88 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
89 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
90 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
91 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
92 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
93 resultMap = Xpetra::toXpetraNonConst(rgTpetraMap);
94#else
95 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Cannot transform Thyra::VectorSpace to Xpetra::Map. This is the general implementation for Tpetra only, but Tpetra is disabled.");
96#endif
97 } // end standard case (no product map)
98 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
99 return resultMap;
100}
101
102template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
103Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
104Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
105 toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
106 using Teuchos::as;
107 using Teuchos::RCP;
108 using Teuchos::rcp_dynamic_cast;
109
110 using ThyMultVecBase = Thyra::MultiVectorBase<Scalar>;
111 using ThyProdMultVecBase = Thyra::ProductMultiVectorBase<Scalar>;
112 using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
113
114 // return value
115 RCP<MultiVector> xpMultVec = Teuchos::null;
116
117 // check whether v is a product multi vector
118 Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase>(v);
119 if (thyProdVec != Teuchos::null) {
120 // SPECIAL CASE: create a nested BlockedMultiVector
121 // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
122 RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
123 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
124
125 // create new Xpetra::BlockedMultiVector
126 xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
127
128 RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec, true);
129
130 // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
131 for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b) {
132 RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
133 // xpBlockMV can be of type MultiVector or BlockedMultiVector
134 RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); // recursive call
135 xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
136 }
137 } else {
138 // STANDARD CASE: no product vector
139#ifdef HAVE_XPETRA_TPETRA
140 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> ConverterT;
141 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpMultVec;
142 typedef Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpTpMultVec;
143 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
144 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
145
146 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
147 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
148 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.");
149 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
150 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
151 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
152 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
153 xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
154#else
155 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Cannot transform Thyra::MultiVector to Xpetra::MultiVector. This is the general implementation for Tpetra only, but Teptra is disabled.");
156#endif
157 } // end standard case
158 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
159 return xpMultVec;
160}
161
162template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
163Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
164Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
165 toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
166 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> cv =
167 Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar>>(v);
168 Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> r =
169 toXpetra(cv, comm);
170 return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(r);
171}
172
173template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
174bool Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
175 isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
176 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(op));
177
178 // check whether we have a Tpetra based Thyra operator
179 bool bIsTpetra = false;
180#ifdef HAVE_XPETRA_TPETRA
181 Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(op);
182 bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
183
184 // for debugging purposes: find out why dynamic cast failed
185 if (!bIsTpetra &&
186#ifdef HAVE_XPETRA_EPETRA
187 Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
188#endif
189 Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar>>(op) == Teuchos::null) {
190 // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
191 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
192 if (Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
193 std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
194 std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
195 std::cout << " properly set!" << std::endl;
196 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
197 }
198 }
199#endif
200
201#if 0
202 // Check whether it is a blocked operator.
203 // If yes, grab the (0,0) block and check the underlying linear algebra
204 // Note: we expect that the (0,0) block exists!
205 if(bIsTpetra == false) {
206 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
207 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
208 if(ThyBlockedOp != Teuchos::null) {
209 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
210 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
211 ThyBlockedOp->getBlock(0,0);
212 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
213 bIsTpetra = isTpetra(b00);
214 }
215 }
216#endif
217
218 return bIsTpetra;
219}
220
221template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
222bool Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
223 isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
224 return false;
225}
226
227template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
228bool Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
229 isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
230 // Check whether it is a blocked operator.
231 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar>> ThyBlockedOp =
232 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar>>(op);
233 if (ThyBlockedOp != Teuchos::null) {
234 return true;
235 }
236 return false;
237}
238
239template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
240Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
241Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
242 toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
243#ifdef HAVE_XPETRA_TPETRA
244 if (isTpetra(op)) {
245 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
246 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
247 // we should also add support for the const versions!
248 // getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
249 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
250 Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
251 Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
252 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat);
253 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
254
255 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
256 Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraNcnstCrsMat));
257 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
258
259 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
260 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
261 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
262 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
263 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
264 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
265 return xpMat;
266 }
267#endif
268
269#ifdef HAVE_XPETRA_EPETRA
270 if (isEpetra(op)) {
271 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
272 }
273#endif
274 return Teuchos::null;
275}
276
277template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
278Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
279Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
280 toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
281 using Teuchos::RCP;
282 using Teuchos::rcp;
283 using Teuchos::rcp_const_cast;
284 using Teuchos::rcp_dynamic_cast;
285
286#ifdef HAVE_XPETRA_TPETRA
287 if (isTpetra(op)) {
288 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
289 typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraOperator_t;
290 RCP<const TpetraOperator_t> TpetraOp = TOE::getConstTpetraOperator(op);
291 TEUCHOS_TEST_FOR_EXCEPT(is_null(TpetraOp));
292 RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat =
293 rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
294 RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat =
295 rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
296
297 RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
298 rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(
299 rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat)));
300 TEUCHOS_TEST_FOR_EXCEPT(is_null(xTpetraCrsMat));
301 RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
302 rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
303 RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
304 rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
305 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
306 rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
307 return xpMat;
308 }
309#endif
310
311#ifdef HAVE_XPETRA_EPETRA
312 if (isEpetra(op)) {
313 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
314 }
315#endif
316 return Teuchos::null;
317}
318
319template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
320Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
321Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
322 toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
323#ifdef HAVE_XPETRA_TPETRA
324 if (isTpetra(op)) {
325 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
326 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
327
328 auto nonConstTpetraOp = Teuchos::rcp_const_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp);
329
330 Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
331 Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(nonConstTpetraOp));
332 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
333
334 Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
335 Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
336 return xpOp;
337 }
338#endif
339
340#ifdef HAVE_XPETRA_EPETRA
341 if (isEpetra(op)) {
342 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
343 }
344#endif
345 return Teuchos::null;
346}
347
348template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
349Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
350Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
351 toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
352#ifdef HAVE_XPETRA_TPETRA
353 if (isTpetra(op)) {
354 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
355 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
356
357 Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
358 Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
359 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
360
361 Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
362 Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
363 return xpOp;
364 }
365#endif
366
367#ifdef HAVE_XPETRA_EPETRA
368 if (isEpetra(op)) {
369 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
370 }
371#endif
372 return Teuchos::null;
373}
374
375template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
376Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
377Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
378 toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op) {
379 using Teuchos::rcp_const_cast;
380 using Teuchos::rcp_dynamic_cast;
381
382 RCP<const Thyra::VectorBase<Scalar>> diag = op->getDiag();
383
384 RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpDiag;
385#ifdef HAVE_XPETRA_TPETRA
386 using thyTpV = Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
387 using tV = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
388 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
389 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
390 if (!tDiag.is_null())
391 xpDiag = Xpetra::toXpetra(tDiag);
392 }
393#endif
394 TEUCHOS_ASSERT(!xpDiag.is_null());
395 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
396 return M;
397}
398
399template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
400Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
401Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
402 toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op) {
403 return toXpetra(Teuchos::rcp_const_cast<Thyra::DiagonalLinearOpBase<Scalar>>(op));
404}
405
406template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
407Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
408Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
409 toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map) {
410 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> thyraMap = Teuchos::null;
411
412 // check whether map is of type BlockedMap
413 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
414 if (bmap.is_null() == false) {
415 Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>> vecSpaces(bmap->getNumMaps());
416 for (size_t i = 0; i < bmap->getNumMaps(); i++) {
417 // we need Thyra GIDs for all the submaps
418 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> vs =
419 Xpetra::ThyraUtils<Scalar, LO, GO, Node>::toThyra(bmap->getMap(i, true));
420 vecSpaces[i] = vs;
421 }
422
423 thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
424 return thyraMap;
425 }
426
427 // standard case
428#ifdef HAVE_XPETRA_TPETRA
429 if (map->lib() == Xpetra::UseTpetra) {
430 Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>> tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(map);
431 if (tpetraMap == Teuchos::null)
432 throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
433 RCP<Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetraMap->getTpetra_Map());
434 thyraMap = thyraTpetraMap;
435 }
436#endif
437
438#ifdef HAVE_XPETRA_EPETRA
439 if (map->lib() == Xpetra::UseEpetra) {
440 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
441 }
442#endif
443
444 return thyraMap;
445}
446
447template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
448Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
449Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
450 toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
451 // create Thyra MultiVector
452#ifdef HAVE_XPETRA_TPETRA
453 if (vec->getMap()->lib() == Xpetra::UseTpetra) {
454 auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
455 RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpMV = Teuchos::rcp_dynamic_cast<const TpetraMultiVector>(vec)->getTpetra_MultiVector();
456 auto thyDomMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(vec->getNumVectors(), vec->getMap()->getComm()));
457 auto thyMV = rcp(new Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
458 thyMV->initialize(thyTpMap, thyDomMap, tpMV);
459 return thyMV;
460 }
461#endif
462
463#ifdef HAVE_XPETRA_EPETRA
464 if (vec->getMap()->lib() == Xpetra::UseEpetra) {
465 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
466 }
467#endif
468
469 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "MultiVector cannot be converted to Thyra.");
470}
471
472template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
473Teuchos::RCP<const Thyra::VectorBase<Scalar>>
474Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
475 toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
476 // create Thyra Vector
477#ifdef HAVE_XPETRA_TPETRA
478 if (vec->getMap()->lib() == Xpetra::UseTpetra) {
479 auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
480 RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpVec = Teuchos::rcp_dynamic_cast<const TpetraVector>(vec)->getTpetra_Vector();
481 auto thyVec = rcp(new Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
482 thyVec->initialize(thyTpMap, tpVec);
483 return thyVec;
484 }
485#endif
486
487#ifdef HAVE_XPETRA_EPETRA
488 if (vec->getMap()->lib() == Xpetra::UseEpetra) {
489 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
490 }
491#endif
492
493 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Vector cannot be converted to Thyra.");
494}
495
496template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
497void Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
498 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) {
499 using Teuchos::as;
500 using Teuchos::RCP;
501 using Teuchos::rcp_dynamic_cast;
502 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
503 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
504 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
505 // typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
506 // typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
507 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
508
509 // copy data from tY_inout to Y_inout
510 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
511 if (prodTarget != Teuchos::null) {
512 RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
513 if (bSourceVec.is_null() == true) {
514 // SPECIAL CASE: target vector is product vector:
515 // update Thyra product multi vector with data from a merged Xpetra multi vector
516
517 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
518 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.");
519
520 for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
521 // access Xpetra data
522 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
523
524 // access Thyra data
525 Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
526 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
527 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
528 const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
529 const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim());
530 RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
531 Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
532
533 // loop over all vectors in multivector
534 for (size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
535 Teuchos::ArrayRCP<const Scalar> xpData = xpSubBlock->getData(j); // access const data from Xpetra object
536
537 // loop over all local rows
538 for (LocalOrdinal i = 0; i < localSubDim; ++i) {
539 (*thyData)(i, j) = xpData[i];
540 }
541 }
542 }
543 } else {
544 // source vector is a blocked multivector
545 // TODO test me
546 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.");
547
548 for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
549 // access Thyra data
550 RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
551
552 Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
553
554 // access Thyra data
555 Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
556 Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
557 }
558 }
559 } else {
560 // STANDARD case:
561 // update Thyra::MultiVector with data from an Xpetra::MultiVector
562
563 // access Thyra data
564 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
565 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
566 const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
567 const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim());
568 RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
569 Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
570
571 // loop over all vectors in multivector
572 for (size_t j = 0; j < source->getNumVectors(); ++j) {
573 Teuchos::ArrayRCP<const Scalar> xpData = source->getData(j); // access const data from Xpetra object
574 // loop over all local rows
575 for (LocalOrdinal i = 0; i < localSubDim; ++i) {
576 (*thyData)(i, j) = xpData[i];
577 }
578 }
579 }
580}
581
582template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
583Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
584Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
585 toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
586 // create a Thyra operator from Xpetra::CrsMatrix
587 Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
588
589 // bool bIsTpetra = false;
590
591#ifdef HAVE_XPETRA_TPETRA
592 Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
593 if (tpetraMat != Teuchos::null) {
594 // bIsTpetra = true;
595 Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
596 Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
597 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
598
599 Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
600 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
601
602 thyraOp = Thyra::createConstLinearOp(tpOperator);
603 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
604 } else {
605#ifdef HAVE_XPETRA_EPETRA
606 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to Tpetra::CrsMatrix failed. Assume matrix should be Epetra then. Epetra needs SC=double, LO=int, and GO=int or GO=long long");
607#else
608 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to Tpetra::CrsMatrix failed. Assume matrix should be Epetra then. No Epetra available");
609#endif
610 }
611 return thyraOp;
612#else
613 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
614 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
615#endif
616}
617
618template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
619Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
620Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
621 toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
622 // create a Thyra operator from Xpetra::CrsMatrix
623 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
624
625 // bool bIsTpetra = false;
626
627#ifdef HAVE_XPETRA_TPETRA
628 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
629 if (tpetraMat != Teuchos::null) {
630 // bIsTpetra = true;
631 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
632 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
633 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
634
635 Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
636 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
637
638 thyraOp = Thyra::createLinearOp(tpOperator);
639 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
640 } else {
641 // cast to TpetraCrsMatrix failed
642#ifdef HAVE_XPETRA_EPETRA
643 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Assuming matrix supposed to be Epetra. Epetra needs SC=double, LO=int, and GO=int or GO=long long");
644#else
645 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Guess, matrix should be Epetra then, but no Epetra available.");
646#endif
647 }
648 return thyraOp;
649#else
650 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
651 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
652#endif
653}
654
655template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
656Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
657Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
658 toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
659 int nRows = mat->Rows();
660 int nCols = mat->Cols();
661
662 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Ablock = mat->getInnermostCrsMatrix();
663 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock);
664 TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
665
666#ifdef HAVE_XPETRA_TPETRA
667 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
668 if (tpetraMat != Teuchos::null) {
669 // create new Thyra blocked operator
670 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar>> blockMat =
671 Thyra::defaultBlockedLinearOp<Scalar>();
672
673 blockMat->beginBlockFill(nRows, nCols);
674
675 for (int r = 0; r < nRows; ++r) {
676 for (int c = 0; c < nCols; ++c) {
677 Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r, c);
678
679 if (xpmat == Teuchos::null) continue; // shortcut for empty blocks
680
681 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thBlock = Teuchos::null;
682
683 // check whether the subblock is again a blocked operator
684 Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpblock =
685 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpmat);
686 if (xpblock != Teuchos::null) {
687 if (xpblock->Rows() == 1 && xpblock->Cols() == 1) {
688 // If it is a single block operator, unwrap it
689 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
690 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
691 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
692 } else {
693 // recursive call for general blocked operators
694 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpblock);
695 }
696 } else {
697 // check whether it is a CRSMatrix object
698 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
699 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
700 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
701 }
702
703 blockMat->setBlock(r, c, thBlock);
704 }
705 }
706
707 blockMat->endBlockFill();
708
709 return blockMat;
710 } else {
711 // tpetraMat == Teuchos::null
712#ifdef HAVE_XPETRA_EPETRA
713 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Assuming matrix supposed to be Epetra. Epetra needs SC=double, LO=int, and GO=int or GO=long long");
714#else
715 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Guess, matrix should be Epetra then, but no Epetra available.");
716#endif
717 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
718 }
719#endif // endif HAVE_XPETRA_TPETRA
720
721 // 4-Aug-2017 JJH Added 2nd condition to avoid "warning: dynamic initialization in unreachable code"
722 // If HAVE_XPETRA_TPETRA is defined, then this method will always return or throw in the if-then-else above.
723#if defined(HAVE_XPETRA_EPETRA) && !defined(HAVE_XPETRA_TPETRA)
724 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
725 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
726#endif // endif HAVE_XPETRA_EPETRA
727}
728
729#ifdef HAVE_XPETRA_EPETRA
730
731#ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
732// implementation of "toThyra" for full specialization on SC=double, LO=GO=int and NO=EpetraNode
733// We need the specialization in the cpp file due to a circle dependency in the .hpp files for BlockedCrsMatrix
734Teuchos::RCP<Thyra::LinearOpBase<double>>
735ThyraUtils<double, int, int, EpetraNode>::toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, int, EpetraNode>>& mat) {
736 int nRows = mat->Rows();
737 int nCols = mat->Cols();
738
739 Teuchos::RCP<Xpetra::Matrix<double, int, int, EpetraNode>> Ablock = mat->getInnermostCrsMatrix();
740 Teuchos::RCP<Xpetra::CrsMatrixWrap<double, int, int, EpetraNode>> Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<double, int, int, EpetraNode>>(Ablock);
741 TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
742
743 bool bTpetra = false;
744 bool bEpetra = false;
745#ifdef HAVE_XPETRA_TPETRA
746 // Note: Epetra is enabled
747#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
748 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
749 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
750 if (tpetraMat != Teuchos::null) bTpetra = true;
751#else
752 bTpetra = false;
753#endif
754#endif
755
756#ifdef HAVE_XPETRA_EPETRA
757 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
758 if (epetraMat != Teuchos::null) bEpetra = true;
759#endif
760
761 TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra
762
763 // create new Thyra blocked operator
764 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar>> blockMat =
765 Thyra::defaultBlockedLinearOp<Scalar>();
766
767 blockMat->beginBlockFill(nRows, nCols);
768
769 for (int r = 0; r < nRows; ++r) {
770 for (int c = 0; c < nCols; ++c) {
771 Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r, c);
772
773 if (xpmat == Teuchos::null) continue; // shortcut for empty blocks
774
775 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thBlock = Teuchos::null;
776
777 // check whether the subblock is again a blocked operator
778 Teuchos::RCP<BlockedCrsMatrix> xpblock = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(xpmat);
779 if (xpblock != Teuchos::null) {
780 if (xpblock->Rows() == 1 && xpblock->Cols() == 1) {
781 // If it is a single block operator, unwrap it
782 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
783 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
784 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
785 } else {
786 // recursive call for general blocked operators
787 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpblock);
788 }
789 } else {
790 // check whether it is a CRSMatrix object
791 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
792 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
793 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
794 }
795
796 blockMat->setBlock(r, c, thBlock);
797 }
798 }
799
800 blockMat->endBlockFill();
801
802 return blockMat;
803}
804#endif // #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
805
806#ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
807// implementation of "toThyra" for full specialization on SC=double, LO=int, GO=long long and NO=EpetraNode
808// We need the specialization in the cpp file due to a circle dependency in the .hpp files for BlockedCrsMatrix
809Teuchos::RCP<Thyra::LinearOpBase<double>>
810ThyraUtils<double, int, long long, EpetraNode>::toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, long long, EpetraNode>>& mat) {
811 int nRows = mat->Rows();
812 int nCols = mat->Cols();
813
814 Teuchos::RCP<Xpetra::Matrix<double, int, long long, EpetraNode>> Ablock = mat->getInnermostCrsMatrix();
815 Teuchos::RCP<Xpetra::CrsMatrixWrap<double, int, long long, EpetraNode>> Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<double, int, long long, EpetraNode>>(Ablock);
816 TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
817
818 bool bTpetra = false;
819 bool bEpetra = false;
820#ifdef HAVE_XPETRA_TPETRA
821 // Note: Epetra is enabled
822#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
823 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
824 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
825 if (tpetraMat != Teuchos::null) bTpetra = true;
826#else
827 bTpetra = false;
828#endif
829#endif
830
831#ifdef HAVE_XPETRA_EPETRA
832 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
833 if (epetraMat != Teuchos::null) bEpetra = true;
834#endif
835
836 TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra
837
838 // create new Thyra blocked operator
839 Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar>> blockMat =
840 Thyra::defaultBlockedLinearOp<Scalar>();
841
842 blockMat->beginBlockFill(nRows, nCols);
843
844 for (int r = 0; r < nRows; ++r) {
845 for (int c = 0; c < nCols; ++c) {
846 Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r, c);
847
848 if (xpmat == Teuchos::null) continue; // shortcut for empty blocks
849
850 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thBlock = Teuchos::null;
851
852 // check whether the subblock is again a blocked operator
853 Teuchos::RCP<BlockedCrsMatrix> xpblock = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(xpmat);
854 if (xpblock != Teuchos::null) {
855 if (xpblock->Rows() == 1 && xpblock->Cols() == 1) {
856 // If it is a single block operator, unwrap it
857 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
858 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
859 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
860 } else {
861 // recursive call for general blocked operators
862 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpblock);
863 }
864 } else {
865 // check whether it is a CRSMatrix object
866 Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
867 TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
868 thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
869 }
870
871 blockMat->setBlock(r, c, thBlock);
872 }
873 }
874
875 blockMat->endBlockFill();
876
877 return blockMat;
878}
879#endif // #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
880
881#endif
882
883} // namespace Xpetra
884
885#endif
886
887#endif
Xpetra_ThyraUtils_decl.hpp
Teuchos::RCP::is_null
bool is_null() const
Teuchos::RCP::ptr
Ptr< T > ptr() const
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition Xpetra_CrsMatrixWrap_decl.hpp:50
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition Xpetra_Exceptions.hpp:65
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::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition Xpetra_MatrixFactory_def.hpp:22
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)
Xpetra::toXpetraNonConst
const RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > toXpetraNonConst(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &map)
Definition Xpetra_TpetraMap_decl.hpp:228
Xpetra::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Definition Xpetra_EpetraCrsGraph.cpp:29