Xpetra Version of the Day
Loading...
Searching...
No Matches
Xpetra_ThyraUtils_decl.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 XPETRA_THYRAUTILS_HPP
11#define XPETRA_THYRAUTILS_HPP
12
13#include "Xpetra_ConfigDefs.hpp"
14#ifdef HAVE_XPETRA_THYRA
15
16#include <typeinfo>
17
18#ifdef HAVE_XPETRA_TPETRA
19#include "Tpetra_ConfigDefs.hpp"
20#endif
21
22#ifdef HAVE_XPETRA_EPETRA
23#include "Epetra_config.h"
24#include "Epetra_CombineMode.h"
25#endif
26
27#include "Xpetra_Map.hpp"
28#include "Xpetra_BlockedMap.hpp"
29#include "Xpetra_BlockedMultiVector.hpp"
30#include "Xpetra_BlockedCrsMatrix.hpp"
31#include "Xpetra_MapUtils.hpp"
32#include "Xpetra_StridedMap.hpp"
33#include "Xpetra_StridedMapFactory.hpp"
34#include "Xpetra_MapExtractor.hpp"
35#include "Xpetra_Matrix.hpp"
36#include "Xpetra_MatrixFactory.hpp"
37#include "Xpetra_CrsMatrixWrap.hpp"
38#include "Xpetra_MultiVectorFactory.hpp"
39
40#include <Thyra_VectorSpaceBase.hpp>
41#include <Thyra_SpmdVectorSpaceBase.hpp>
42#include <Thyra_ProductVectorSpaceBase.hpp>
43#include <Thyra_ProductMultiVectorBase.hpp>
44#include <Thyra_VectorSpaceBase.hpp>
45#include <Thyra_DefaultProductVectorSpace.hpp>
46#include <Thyra_DefaultBlockedLinearOp.hpp>
47#include <Thyra_LinearOpBase.hpp>
48#include "Thyra_DiagonalLinearOpBase.hpp"
49#include <Thyra_DetachedMultiVectorView.hpp>
50#include <Thyra_MultiVectorStdOps.hpp>
51
52#ifdef HAVE_XPETRA_TPETRA
53#include <Thyra_TpetraThyraWrappers.hpp>
54#include <Thyra_TpetraVector.hpp>
55#include <Thyra_TpetraMultiVector.hpp>
56#include <Thyra_TpetraVectorSpace.hpp>
57#include <Tpetra_Map.hpp>
58#include <Tpetra_Vector.hpp>
59#include <Tpetra_CrsMatrix.hpp>
60#include <Xpetra_TpetraMap.hpp>
61#include <Xpetra_TpetraOperator.hpp>
62#include <Xpetra_TpetraCrsMatrix.hpp>
63#endif
64#ifdef HAVE_XPETRA_EPETRA
65#include <Thyra_EpetraLinearOp.hpp>
66#include <Thyra_EpetraThyraWrappers.hpp>
67#include <Thyra_SpmdVectorBase.hpp>
68#include <Thyra_get_Epetra_Operator.hpp>
69#include <Epetra_Map.h>
70#include <Epetra_Vector.h>
71#include <Epetra_CrsMatrix.h>
72#include <Xpetra_EpetraMap.hpp>
73#include <Xpetra_EpetraCrsMatrix.hpp>
74#endif
75
76namespace Xpetra {
77
78template <class Scalar,
79 class LocalOrdinal = int,
80 class GlobalOrdinal = LocalOrdinal,
81 class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
82class ThyraUtils {
83 private:
84#undef XPETRA_THYRAUTILS_SHORT
85#include "Xpetra_UseShortNames.hpp"
86
87 public:
88 static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>
89 toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm, std::vector<size_t>& stridingInfo, LocalOrdinal stridedBlockId = -1, GlobalOrdinal offset = 0);
90
91 static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
92 toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm);
93
94 // const version
95 static Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
96 toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm);
97
98 // non-const version
99 static Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
100 toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm);
101
102 static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op);
103
104 static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op);
105
106 static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op);
107
108 static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
109 toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op);
110
111 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
112 toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op);
113
114 static Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
115 toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op);
116
117 static Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
118 toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op);
119
120 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
121 toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op);
122
123 static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
124 toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op);
125
126 static Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
127 toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map);
128
129 static Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
130 toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec);
131
132 static Teuchos::RCP<const Thyra::VectorBase<Scalar>>
133 toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec);
134
135 // update Thyra multi vector with data from Xpetra multi vector
136 // In case of a Thyra::ProductMultiVector the Xpetra::MultiVector is splitted into its subparts using a provided MapExtractor
137 static void
138 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);
139
140 static Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
141 toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat);
142
143 static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
144 toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat);
145
146 static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
147 toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat);
148
149}; // end Utils class
150
151// full specialization for Epetra support
152// Note, that Thyra only has support for Epetra (not for Epetra64)
153#ifdef HAVE_XPETRA_EPETRA
154
155#ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
156template <>
157class ThyraUtils<double, int, int, EpetraNode> {
158 public:
159 typedef double Scalar;
160 typedef int LocalOrdinal;
161 typedef int GlobalOrdinal;
162 typedef EpetraNode Node;
163
164 private:
165#undef XPETRA_THYRAUTILS_SHORT
166#include "Xpetra_UseShortNames.hpp"
167
168 public:
169 static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>
170 toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm, std::vector<size_t>& stridingInfo, LocalOrdinal stridedBlockId = -1, GlobalOrdinal offset = 0) {
171 Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map = ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(vectorSpace, comm);
172
173 if (stridedBlockId == -1) {
174 TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo.size() != 0);
175 } else {
176 TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo[stridedBlockId] != 0);
177 }
178
179 Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>> ret = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(map, stridingInfo, stridedBlockId, offset);
180 return ret;
181 }
182
183 static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
184 toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
185 using Teuchos::as;
186 using Teuchos::RCP;
187 using Teuchos::rcp_dynamic_cast;
188 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
189 typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
190 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyUtils;
191
192 RCP<Map> resultMap = Teuchos::null;
193
194 RCP<const ThyProdVecSpaceBase> prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase>(vectorSpace);
195 if (prodVectorSpace != Teuchos::null) {
196 // SPECIAL CASE: product Vector space
197 // collect all submaps to store them in a hierarchical BlockedMap object
198 TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks() == 0, std::logic_error, "Found a product vector space with zero blocks.");
199 std::vector<RCP<const Map>> mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
200 std::vector<RCP<const Map>> mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
201 for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
202 RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
203 // can be of type Map or BlockedMap (containing Thyra GIDs)
204 mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
205 }
206
207 // get offsets for submap GIDs
208 // we need that for the full map (Xpetra GIDs)
209 std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(), 0);
210 for (int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
211 gidOffsets[i] = mapsThyra[i - 1]->getMaxAllGlobalIndex() + gidOffsets[i - 1] + 1;
212 }
213
214 for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
215 RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
216 // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
217 mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
218 }
219
220 resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapsXpetra, mapsThyra));
221 } else {
222 // STANDARD CASE: no product map
223 // Epetra/Tpetra specific code to access the underlying map data
224
225 // check whether we have a Tpetra based Thyra operator
226 bool bIsTpetra = false;
227#ifdef HAVE_XPETRA_TPETRA
228#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
229 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
230 Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vectorSpace);
231 bIsTpetra = Teuchos::is_null(tpetra_vsc) ? false : true;
232#endif
233#endif
234
235 // check whether we have an Epetra based Thyra operator
236 bool bIsEpetra = !bIsTpetra; // note: this is a little bit fragile!
237
238#ifdef HAVE_XPETRA_TPETRA
239 if (bIsTpetra) {
240#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
241 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
242 typedef Thyra::VectorBase<Scalar> ThyVecBase;
243 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> TpMap;
244 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpVector;
245 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
246 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
247 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
248 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
249 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
250 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
251 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
252
253 resultMap = Xpetra::toXpetraNonConst(rgTpetraMap);
254#else
255 throw Xpetra::Exceptions::RuntimeError("Problem AAA. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
256#endif
257 }
258#endif
259
260#ifdef HAVE_XPETRA_EPETRA
261 if (bIsEpetra) {
262 // RCP<const Epetra_Map> epMap = Teuchos::null;
263 RCP<const Epetra_Map>
264 epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map>>(vectorSpace, "epetra_map");
265 if (!Teuchos::is_null(epetra_map)) {
266 resultMap = Teuchos::rcp(new Xpetra::EpetraMapT<GlobalOrdinal, Node>(epetra_map));
267 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
268 } else {
269 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "No Epetra_Map data found in Thyra::VectorSpace.");
270 }
271 }
272#endif
273 } // end standard case (no product map)
274 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
275 return resultMap;
276 }
277
278 // const version
279 static Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
280 toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
281 using Teuchos::as;
282 using Teuchos::RCP;
283 using Teuchos::rcp_dynamic_cast;
284
285 using ThyProdMultVecBase = Thyra::ProductMultiVectorBase<Scalar>;
286 using ThyMultVecBase = Thyra::MultiVectorBase<Scalar>;
287 using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
288
289 // return value
290 RCP<MultiVector> xpMultVec = Teuchos::null;
291
292 // check whether v is a product multi vector
293 Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase>(v);
294 if (thyProdVec != Teuchos::null) {
295 // SPECIAL CASE: create a nested BlockedMultiVector
296 // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
297 RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
298 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
299
300 // create new Xpetra::BlockedMultiVector
301 xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
302
303 RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec, true);
304
305 // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
306 for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b) {
307 RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
308 // xpBlockMV can be of type MultiVector or BlockedMultiVector
309 RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); // recursive call
310 xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
311 }
312
313 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
314 return xpMultVec;
315 } else {
316 // STANDARD CASE: no product vector
317 // Epetra/Tpetra specific code to access the underlying map data
318 bool bIsTpetra = false;
319#ifdef HAVE_XPETRA_TPETRA
320#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
321 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
322
323 // typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
324 // typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
325 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
326 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> ConverterT;
327 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpMultVec;
328 typedef Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpTpMultVec;
329 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
330
331 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
332 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
333 if (thyraTpetraMultiVector != Teuchos::null) {
334 bIsTpetra = true;
335 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
336 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
337 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
338 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
339 xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
340 }
341#endif
342#endif
343
344#ifdef HAVE_XPETRA_EPETRA
345 if (bIsTpetra == false) {
346 // no product vector
347 Teuchos::RCP<Map> map = ThyUtils::toXpetra(v->range(), comm);
348 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(map));
349 RCP<Xpetra::EpetraMapT<GlobalOrdinal, Node>> xeMap = rcp_dynamic_cast<Xpetra::EpetraMapT<GlobalOrdinal, Node>>(map, true);
350 RCP<const Epetra_Map> eMap = xeMap->getEpetra_MapRCP();
351 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(eMap));
352 Teuchos::RCP<const Epetra_MultiVector> epMultVec = Thyra::get_Epetra_MultiVector(*eMap, v);
353 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epMultVec));
354 RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<Epetra_MultiVector>(epMultVec);
355 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epNonConstMultVec));
356 xpMultVec = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>(epNonConstMultVec));
357 }
358#endif
359 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
360 return xpMultVec;
361 } // end standard case
362 }
363
364 // non-const version
365 static Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
366 toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
367 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> cv =
368 Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar>>(v);
369 Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> r =
370 toXpetra(cv, comm);
371 return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(r);
372 }
373
374 static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
375 // check whether we have a Tpetra based Thyra operator
376 bool bIsTpetra = false;
377#ifdef HAVE_XPETRA_TPETRA
378#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
379 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
380
381 Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(op);
382 bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
383
384 // for debugging purposes: find out why dynamic cast failed
385 if (!bIsTpetra &&
386#ifdef HAVE_XPETRA_EPETRA
387 Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
388#endif
389 Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar>>(op) == Teuchos::null) {
390 // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
391 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
392 if (Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
393 std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
394 std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
395 std::cout << " properly set!" << std::endl;
396 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
397 }
398 }
399#endif
400#endif
401
402#if 0
403 // Check whether it is a blocked operator.
404 // If yes, grab the (0,0) block and check the underlying linear algebra
405 // Note: we expect that the (0,0) block exists!
406 if(bIsTpetra == false) {
407 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
408 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
409 if(ThyBlockedOp != Teuchos::null) {
410 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
411 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
412 ThyBlockedOp->getBlock(0,0);
413 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
414 bIsTpetra = isTpetra(b00);
415 }
416 }
417#endif
418
419 return bIsTpetra;
420 }
421
422 static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
423 // check whether we have an Epetra based Thyra operator
424 bool bIsEpetra = false;
425
426#ifdef HAVE_XPETRA_EPETRA
427 Teuchos::RCP<const Thyra::EpetraLinearOp> epetra_op = Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op, false);
428 bIsEpetra = Teuchos::is_null(epetra_op) ? false : true;
429#endif
430
431#if 0
432 // Check whether it is a blocked operator.
433 // If yes, grab the (0,0) block and check the underlying linear algebra
434 // Note: we expect that the (0,0) block exists!
435 if(bIsEpetra == false) {
436 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
437 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op,false);
438 if(ThyBlockedOp != Teuchos::null) {
439 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
440 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
441 ThyBlockedOp->getBlock(0,0);
442 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
443 bIsEpetra = isEpetra(b00);
444 }
445 }
446#endif
447
448 return bIsEpetra;
449 }
450
451 static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
452 // Check whether it is a blocked operator.
453 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar>> ThyBlockedOp =
454 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar>>(op);
455 if (ThyBlockedOp != Teuchos::null) {
456 return true;
457 }
458 return false;
459 }
460
461 static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
462 toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
463#ifdef HAVE_XPETRA_TPETRA
464 if (isTpetra(op)) {
465#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
466 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
467
468 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
469 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
470 // we should also add support for the const versions!
471 // getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
472 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
473 Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp);
474 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraRowMat));
475 Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
476 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat);
477 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
478
479 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
480 Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraNcnstCrsMat));
481 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
482 Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> ret =
483 Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat);
484 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
485 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
486 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
487 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
488 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
489 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
490 return xpMat;
491#else
492 throw Xpetra::Exceptions::RuntimeError("Problem BBB. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
493#endif
494 }
495#endif
496
497#ifdef HAVE_XPETRA_EPETRA
498 if (isEpetra(op)) {
499 Teuchos::RCP<const Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*op);
500 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
501 Teuchos::RCP<const Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_op, true);
502 Teuchos::RCP<const Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(epetra_rowmat, true);
503 Teuchos::RCP<Epetra_CrsMatrix> epetra_ncnstcrsmat = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(epetra_crsmat);
504 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_ncnstcrsmat));
505
506 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpetraCrsMat =
507 Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>(epetra_ncnstcrsmat));
508 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
509
510 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
511 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xEpetraCrsMat, true);
512 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
513 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
514 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
515 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
516 return xpMat;
517 }
518#endif
519 return Teuchos::null;
520 }
521
522 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
523 toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
524#ifdef HAVE_XPETRA_TPETRA
525 if (isTpetra(op)) {
526#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
527 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
528
529 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
530 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
531 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
532 Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
533 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
534
535 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
536 Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraCrsMat));
537 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
538 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
539 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
540 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
541 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
542 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
543 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
544 return xpMat;
545#else
546 throw Xpetra::Exceptions::RuntimeError("Problem CCC. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
547#endif
548 }
549#endif
550
551#ifdef HAVE_XPETRA_EPETRA
552 if (isEpetra(op)) {
553 Teuchos::RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*op);
554 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
555 Teuchos::RCP<Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetra_op, true);
556 Teuchos::RCP<Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetra_rowmat, true);
557
558 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpetraCrsMat =
559 Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>(epetra_crsmat));
560 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
561 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
562 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xEpetraCrsMat, true);
563 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
564 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
565 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
566 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
567 return xpMat;
568 }
569#endif
570 return Teuchos::null;
571 }
572
573 static Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
574 toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
575 return toXpetraOperator(Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar>>(op));
576
577 // #ifdef HAVE_XPETRA_TPETRA
578 // if(isTpetra(op)) {
579 // typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
580 // Teuchos::RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getConstTpetraOperator(op);
581
582 // Teuchos::RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > nonConstTpetraOp =
583 // Teuchos::rcp_const_cast<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
584
585 // Teuchos::RCP<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xTpetraOp =
586 // Teuchos::rcp(new Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>(nonConstTpetraOp));
587 // TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
588
589 // Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xpOp =
590 // Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xTpetraOp, true);
591 // return xpOp;
592 // }
593 // #endif
594
595 // #ifdef HAVE_XPETRA_EPETRA
596 // if(isEpetra(op)) {
597 // TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
598 // }
599 // #endif
600 // return Teuchos::null;
601 }
602
603 static Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
604 toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
605#ifdef HAVE_XPETRA_TPETRA
606 if (isTpetra(op)) {
607 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
608 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
609
610 Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
611 Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
612 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
613
614 Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
615 Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
616 return xpOp;
617 }
618#endif
619
620#ifdef HAVE_XPETRA_EPETRA
621 if (isEpetra(op)) {
622 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
623 }
624#endif
625 return Teuchos::null;
626 }
627
628 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
629 toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op) {
630 using Teuchos::rcp_const_cast;
631 using Teuchos::rcp_dynamic_cast;
632
633 RCP<const Thyra::VectorBase<Scalar>> diag = op->getDiag();
634
635 RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpDiag;
636#ifdef HAVE_XPETRA_TPETRA
637 using thyTpV = Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
638 using tV = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
639 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
640 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
641 if (!tDiag.is_null())
642 xpDiag = Xpetra::toXpetra(tDiag);
643 }
644#endif
645#ifdef HAVE_XPETRA_EPETRA
646 using ThyVSBase = Thyra::SpmdVectorSpaceBase<Scalar>;
647 if (xpDiag.is_null()) {
648 RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(op->range())->getComm());
649 RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(op->range()), comm);
650 if (!map.is_null()) {
651 RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
652 RCP<Epetra_Vector> nceDiag = rcp_const_cast<Epetra_Vector>(eDiag);
653 RCP<Xpetra::EpetraVectorT<int, Node>> xpEpDiag = rcp(new Xpetra::EpetraVectorT<int, Node>(nceDiag));
654 xpDiag = rcp_dynamic_cast<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpEpDiag, true);
655 }
656 }
657#endif
658 TEUCHOS_ASSERT(!xpDiag.is_null());
659 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
660 return M;
661 }
662
663 static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
664 toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op) {
665 return toXpetra(Teuchos::rcp_const_cast<Thyra::DiagonalLinearOpBase<Scalar>>(op));
666 }
667
668 static Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
669 toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map) {
670 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> thyraMap = Teuchos::null;
671
672 // check whether map is of type BlockedMap
673 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
674 if (bmap.is_null() == false) {
675 Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>> vecSpaces(bmap->getNumMaps());
676 for (size_t i = 0; i < bmap->getNumMaps(); i++) {
677 // we need Thyra GIDs for all the submaps
678 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> vs =
679 Xpetra::ThyraUtils<Scalar, LO, GO, Node>::toThyra(bmap->getMap(i, true));
680 vecSpaces[i] = vs;
681 }
682
683 thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
684 return thyraMap;
685 }
686
687 // standard case
688#ifdef HAVE_XPETRA_TPETRA
689 if (map->lib() == Xpetra::UseTpetra) {
690#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
691 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
692 Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>> tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(map);
693 if (tpetraMap == Teuchos::null)
694 throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
695 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tpMap = tpetraMap->getTpetra_Map();
696 RCP<Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
697 thyraMap = thyraTpetraMap;
698#else
699 throw Xpetra::Exceptions::RuntimeError("Problem DDD. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
700#endif
701 }
702#endif
703
704#ifdef HAVE_XPETRA_EPETRA
705 if (map->lib() == Xpetra::UseEpetra) {
706 Teuchos::RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>> epetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(map);
707 if (epetraMap == Teuchos::null)
708 throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
709 RCP<const Thyra::VectorSpaceBase<Scalar>> thyraEpetraMap = Thyra::create_VectorSpace(epetraMap->getEpetra_MapRCP());
710 thyraMap = thyraEpetraMap;
711 }
712#endif
713
714 return thyraMap;
715 }
716
717 static Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
718 toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
719 // create Thyra MultiVector
720#ifdef HAVE_XPETRA_TPETRA
721 if (vec->getMap()->lib() == Xpetra::UseTpetra) {
722 auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
723 RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpMV = Teuchos::rcp_dynamic_cast<const TpetraMultiVector>(vec)->getTpetra_MultiVector();
724 auto thyDomMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(vec->getNumVectors(), vec->getMap()->getComm()));
725 auto thyMV = rcp(new Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
726 thyMV->initialize(thyTpMap, thyDomMap, tpMV);
727 return thyMV;
728 }
729#endif
730
731#ifdef HAVE_XPETRA_EPETRA
732 if (vec->getMap()->lib() == Xpetra::UseEpetra) {
733 auto thyEpMap = Thyra::create_VectorSpace(Teuchos::rcp_dynamic_cast<const EpetraMapT<GlobalOrdinal, Node>>(vec->getMap())->getEpetra_MapRCP());
734 auto epMV = Teuchos::rcp_dynamic_cast<const EpetraMultiVectorT<GlobalOrdinal, Node>>(vec)->getEpetra_MultiVector();
735 auto thyMV = Thyra::create_MultiVector(epMV, thyEpMap);
736 return thyMV;
737 }
738#endif
739
740 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "MultiVector cannot be converted to Thyra.");
741 }
742
743 static Teuchos::RCP<const Thyra::VectorBase<Scalar>>
744 toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
745 // create Thyra Vector
746#ifdef HAVE_XPETRA_TPETRA
747 if (vec->getMap()->lib() == Xpetra::UseTpetra) {
748 auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
749 RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpVec = Teuchos::rcp_dynamic_cast<const TpetraVector>(vec)->getTpetra_Vector();
750 auto thyVec = rcp(new Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
751 thyVec->initialize(thyTpMap, tpVec);
752 return thyVec;
753 }
754#endif
755
756#ifdef HAVE_XPETRA_EPETRA
757 if (vec->getMap()->lib() == Xpetra::UseEpetra) {
758 auto thyEpMap = Thyra::create_VectorSpace(Teuchos::rcp_dynamic_cast<const EpetraMapT<GlobalOrdinal, Node>>(vec->getMap())->getEpetra_MapRCP());
759 auto epVec = rcp(Teuchos::rcp_dynamic_cast<const EpetraVectorT<GlobalOrdinal, Node>>(vec)->getEpetra_Vector(), false);
760 auto thyVec = Thyra::create_Vector(epVec, thyEpMap);
761 return thyVec;
762 }
763#endif
764
765 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Vector cannot be converted to Thyra.");
766 }
767
768 static void 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) {
769 using Teuchos::as;
770 using Teuchos::RCP;
771 using Teuchos::rcp_dynamic_cast;
772 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
773 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
774 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
775 // typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
776 // typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
777 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
778
779 // copy data from tY_inout to Y_inout
780 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
781 if (prodTarget != Teuchos::null) {
782 RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
783 if (bSourceVec.is_null() == true) {
784 // SPECIAL CASE: target vector is product vector:
785 // update Thyra product multi vector with data from a merged Xpetra multi vector
786
787 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
788 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.");
789
790 for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
791 // access Xpetra data
792 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
793
794 // access Thyra data
795 Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
796 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
797 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
798 const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
799 const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim());
800 RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
801 Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
802
803 // loop over all vectors in multivector
804 for (size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
805 Teuchos::ArrayRCP<const Scalar> xpData = xpSubBlock->getData(j); // access const data from Xpetra object
806
807 // loop over all local rows
808 for (LocalOrdinal i = 0; i < localSubDim; ++i) {
809 (*thyData)(i, j) = xpData[i];
810 }
811 }
812 }
813 } else {
814 // source vector is a blocked multivector
815 // TODO test me
816 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.");
817
818 for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
819 // access Thyra data
820 RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
821
822 Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
823
824 // access Thyra data
825 Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
826 Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
827 }
828 }
829 } else {
830 // STANDARD case:
831 // update Thyra::MultiVector with data from an Xpetra::MultiVector
832
833 // access Thyra data
834 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
835 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
836 const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
837 const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim());
838 RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
839 Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
840
841 // loop over all vectors in multivector
842 for (size_t j = 0; j < source->getNumVectors(); ++j) {
843 Teuchos::ArrayRCP<const Scalar> xpData = source->getData(j); // access const data from Xpetra object
844 // loop over all local rows
845 for (LocalOrdinal i = 0; i < localSubDim; ++i) {
846 (*thyData)(i, j) = xpData[i];
847 }
848 }
849 }
850 }
851
852 static Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
853 toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
854 // create a Thyra operator from Xpetra::CrsMatrix
855 Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
856
857#ifdef HAVE_XPETRA_TPETRA
858 Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
859 if (tpetraMat != Teuchos::null) {
860#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
861 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
862
863 Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
864 Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
865 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
866
867 Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
868 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
869
870 thyraOp = Thyra::createConstLinearOp(tpOperator);
871 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
872#else
873 throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
874#endif
875 }
876#endif
877
878#ifdef HAVE_XPETRA_EPETRA
879 Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
880 if (epetraMat != Teuchos::null) {
881 Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
882 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpCrsMat));
883 Teuchos::RCP<const Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrix();
884 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
885
886 Teuchos::RCP<const Thyra::EpetraLinearOp> thyraEpOp = Thyra::epetraLinearOp(epCrsMat, "op");
887 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
888 thyraOp = thyraEpOp;
889 }
890#endif
891 return thyraOp;
892 }
893
894 static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
895 toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
896 // create a Thyra operator from Xpetra::CrsMatrix
897 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
898
899#ifdef HAVE_XPETRA_TPETRA
900 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
901 if (tpetraMat != Teuchos::null) {
902#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
903 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
904
905 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
906 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
907 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
908
909 Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
910 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
911
912 thyraOp = Thyra::createLinearOp(tpOperator);
913 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
914#else
915 throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
916#endif
917 }
918#endif
919
920#ifdef HAVE_XPETRA_EPETRA
921 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
922 if (epetraMat != Teuchos::null) {
923 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat, true);
924 Teuchos::RCP<Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrixNonConst();
925 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
926
927 Teuchos::RCP<Thyra::EpetraLinearOp> thyraEpOp = Thyra::nonconstEpetraLinearOp(epCrsMat, "op");
928 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
929 thyraOp = thyraEpOp;
930 }
931#endif
932 return thyraOp;
933 }
934
935 static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
936 toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, int, EpetraNode>>& mat);
937
938}; // specialization on SC=double, LO=GO=int and NO=EpetraNode
939#endif // #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
940
941#ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
942template <>
943class ThyraUtils<double, int, long long, EpetraNode> {
944 public:
945 typedef double Scalar;
946 typedef int LocalOrdinal;
947 typedef long long GlobalOrdinal;
948 typedef EpetraNode Node;
949
950 private:
951#undef XPETRA_THYRAUTILS_SHORT
952#include "Xpetra_UseShortNames.hpp"
953
954 public:
955 static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>
956 toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm, std::vector<size_t>& stridingInfo, LocalOrdinal stridedBlockId = -1, GlobalOrdinal offset = 0) {
957 Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map = ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(vectorSpace, comm);
958
959 if (stridedBlockId == -1) {
960 TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo.size() != 0);
961 } else {
962 TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo[stridedBlockId] != 0);
963 }
964
965 Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>> ret = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(map, stridingInfo, stridedBlockId, offset);
966 return ret;
967 }
968
969 static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
970 toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
971 using Teuchos::as;
972 using Teuchos::RCP;
973 using Teuchos::rcp_dynamic_cast;
974 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
975 typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
976 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyUtils;
977
978 RCP<const ThyProdVecSpaceBase> prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase>(vectorSpace);
979 if (prodVectorSpace != Teuchos::null) {
980 // SPECIAL CASE: product Vector space
981 // collect all submaps to store them in a hierarchical BlockedMap object
982 TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks() == 0, std::logic_error, "Found a product vector space with zero blocks.");
983 std::vector<RCP<const Map>> mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
984 std::vector<RCP<const Map>> mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
985 for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
986 RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
987 // can be of type Map or BlockedMap (containing Thyra GIDs)
988 mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
989 }
990
991 // get offsets for submap GIDs
992 // we need that for the full map (Xpetra GIDs)
993 std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(), 0);
994 for (int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
995 gidOffsets[i] = mapsThyra[i - 1]->getMaxAllGlobalIndex() + gidOffsets[i - 1] + 1;
996 }
997
998 for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
999 RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
1000 // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
1001 mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
1002 }
1003
1004 Teuchos::RCP<Map> resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapsXpetra, mapsThyra));
1005 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
1006 return resultMap;
1007 } else {
1008 // STANDARD CASE: no product map
1009 // Epetra/Tpetra specific code to access the underlying map data
1010
1011 // check whether we have a Tpetra based Thyra operator
1012 bool bIsTpetra = false;
1013#ifdef HAVE_XPETRA_TPETRA
1014#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1015 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1016 Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vectorSpace);
1017 bIsTpetra = Teuchos::is_null(tpetra_vsc) ? false : true;
1018#endif
1019#endif
1020
1021 // check whether we have an Epetra based Thyra operator
1022 bool bIsEpetra = !bIsTpetra; // note: this is a little bit fragile!
1023
1024#ifdef HAVE_XPETRA_TPETRA
1025 if (bIsTpetra) {
1026#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1027 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1028 typedef Thyra::VectorBase<Scalar> ThyVecBase;
1029 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> TpMap;
1030 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpVector;
1031 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1032 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
1033 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
1034 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
1035 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
1036 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
1037 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
1038
1039 RCP<Map> rgXpetraMap = Xpetra::toXpetraNonConst(rgTpetraMap);
1040 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgXpetraMap));
1041 return rgXpetraMap;
1042#else
1043 throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1044#endif
1045 }
1046#endif
1047
1048#ifdef HAVE_XPETRA_EPETRA
1049 if (bIsEpetra) {
1050 // RCP<const Epetra_Map> epMap = Teuchos::null;
1051 RCP<const Epetra_Map>
1052 epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map>>(vectorSpace, "epetra_map");
1053 if (!Teuchos::is_null(epetra_map)) {
1054 Teuchos::RCP<Map> rgXpetraMap = Teuchos::rcp(new Xpetra::EpetraMapT<GlobalOrdinal, Node>(epetra_map));
1055 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgXpetraMap));
1056 return rgXpetraMap;
1057 } else {
1058 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "No Epetra_Map data found in Thyra::VectorSpace.");
1059 }
1060 }
1061#endif
1062 } // end standard case (no product map)
1063 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Cannot transform Thyra::VectorSpace to Xpetra::Map.");
1064 // return Teuchos::null; // unreachable
1065 }
1066
1067 // const version
1068 static Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1069 toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
1070 using Teuchos::as;
1071 using Teuchos::RCP;
1072 using Teuchos::rcp_dynamic_cast;
1073 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1074 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1075 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyUtils;
1076
1077 // return value
1078 RCP<MultiVector> xpMultVec = Teuchos::null;
1079
1080 // check whether v is a product multi vector
1081 Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase>(v);
1082 if (thyProdVec != Teuchos::null) {
1083 // SPECIAL CASE: create a nested BlockedMultiVector
1084 // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
1085 RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
1086 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
1087
1088 // create new Xpetra::BlockedMultiVector
1089 xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
1090
1091 RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec);
1092 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpBlockedMultVec));
1093
1094 // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
1095 for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b) {
1096 RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
1097 // xpBlockMV can be of type MultiVector or BlockedMultiVector
1098 RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); // recursive call
1099 xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
1100 }
1101
1102 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1103 return xpMultVec;
1104 } else {
1105 // STANDARD CASE: no product vector
1106 // Epetra/Tpetra specific code to access the underlying map data
1107 bool bIsTpetra = false;
1108#ifdef HAVE_XPETRA_TPETRA
1109#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1110 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1111
1112 // typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
1113 // typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
1114 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
1115 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> ConverterT;
1116 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpMultVec;
1117 typedef Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpTpMultVec;
1118 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
1119
1120 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
1121 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
1122 if (thyraTpetraMultiVector != Teuchos::null) {
1123 bIsTpetra = true;
1124 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
1125 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
1126 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
1127 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
1128 xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
1129 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1130 return xpMultVec;
1131 }
1132#endif
1133#endif
1134
1135#ifdef HAVE_XPETRA_EPETRA
1136 if (bIsTpetra == false) {
1137 // no product vector
1138 Teuchos::RCP<Map> map = ThyUtils::toXpetra(v->range(), comm);
1139 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(map));
1140 RCP<Xpetra::EpetraMapT<GlobalOrdinal, Node>> xeMap = rcp_dynamic_cast<Xpetra::EpetraMapT<GlobalOrdinal, Node>>(map, true);
1141 RCP<const Epetra_Map> eMap = xeMap->getEpetra_MapRCP();
1142 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(eMap));
1143 Teuchos::RCP<const Epetra_MultiVector> epMultVec = Thyra::get_Epetra_MultiVector(*eMap, v);
1144 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epMultVec));
1145 RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<Epetra_MultiVector>(epMultVec);
1146 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epNonConstMultVec));
1147 xpMultVec = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>(epNonConstMultVec));
1148 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1149 return xpMultVec;
1150 }
1151#endif
1152 } // end standard case
1153 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Cannot transform Thyra::MultiVector to Xpetra::MultiVector.");
1154 // return Teuchos::null; // unreachable
1155 }
1156
1157 // non-const version
1158 static Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1159 toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
1160 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> cv =
1161 Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar>>(v);
1162 Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> r =
1163 toXpetra(cv, comm);
1164 return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(r);
1165 }
1166
1167 static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1168 // check whether we have a Tpetra based Thyra operator
1169 bool bIsTpetra = false;
1170#ifdef HAVE_XPETRA_TPETRA
1171#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1172 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1173
1174 Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(op);
1175 bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
1176
1177 // for debugging purposes: find out why dynamic cast failed
1178 if (!bIsTpetra &&
1179#ifdef HAVE_XPETRA_EPETRA
1180 Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
1181#endif
1182 Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar>>(op) == Teuchos::null) {
1183 // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
1184 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
1185 if (Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
1186 std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
1187 std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
1188 std::cout << " properly set!" << std::endl;
1189 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
1190 }
1191 }
1192#endif
1193#endif
1194
1195#if 0
1196 // Check whether it is a blocked operator.
1197 // If yes, grab the (0,0) block and check the underlying linear algebra
1198 // Note: we expect that the (0,0) block exists!
1199 if(bIsTpetra == false) {
1200 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1201 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
1202 if(ThyBlockedOp != Teuchos::null) {
1203 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1204 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1205 ThyBlockedOp->getBlock(0,0);
1206 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1207 bIsTpetra = isTpetra(b00);
1208 }
1209 }
1210#endif
1211
1212 return bIsTpetra;
1213 }
1214
1215 static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1216 // check whether we have an Epetra based Thyra operator
1217 bool bIsEpetra = false;
1218
1219#ifdef HAVE_XPETRA_EPETRA
1220 Teuchos::RCP<const Thyra::EpetraLinearOp> epetra_op = Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op, false);
1221 bIsEpetra = Teuchos::is_null(epetra_op) ? false : true;
1222#endif
1223
1224#if 0
1225 // Check whether it is a blocked operator.
1226 // If yes, grab the (0,0) block and check the underlying linear algebra
1227 // Note: we expect that the (0,0) block exists!
1228 if(bIsEpetra == false) {
1229 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1230 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op,false);
1231 if(ThyBlockedOp != Teuchos::null) {
1232 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1233 Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1234 ThyBlockedOp->getBlock(0,0);
1235 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1236 bIsEpetra = isEpetra(b00);
1237 }
1238 }
1239#endif
1240
1241 return bIsEpetra;
1242 }
1243
1244 static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1245 // Check whether it is a blocked operator.
1246 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar>> ThyBlockedOp =
1247 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar>>(op);
1248 if (ThyBlockedOp != Teuchos::null) {
1249 return true;
1250 }
1251 return false;
1252 }
1253
1254 static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1255 toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1256#ifdef HAVE_XPETRA_TPETRA
1257 if (isTpetra(op)) {
1258#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1259 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1260
1261 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1262 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
1263 // we should also add support for the const versions!
1264 // getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
1265 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1266 Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
1267 Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
1268 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat);
1269 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
1270
1271 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
1272 Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraNcnstCrsMat));
1273 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1274
1275 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1276 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
1277 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1278 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1279 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1280 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1281 return xpMat;
1282#else
1283 throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1284#endif
1285 }
1286#endif
1287
1288#ifdef HAVE_XPETRA_EPETRA
1289 if (isEpetra(op)) {
1290 Teuchos::RCP<const Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*op);
1291 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1292 Teuchos::RCP<const Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_op, true);
1293 Teuchos::RCP<const Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(epetra_rowmat, true);
1294 Teuchos::RCP<Epetra_CrsMatrix> epetra_ncnstcrsmat = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(epetra_crsmat);
1295 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_ncnstcrsmat));
1296
1297 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpetraCrsMat =
1298 Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>(epetra_ncnstcrsmat));
1299 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1300
1301 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1302 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xEpetraCrsMat, true);
1303 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1304 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1305 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1306 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1307 return xpMat;
1308 }
1309#endif
1310 return Teuchos::null;
1311 }
1312
1313 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1314 toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
1315#ifdef HAVE_XPETRA_TPETRA
1316 if (isTpetra(op)) {
1317#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1318 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1319
1320 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1321 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
1322 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1323 Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
1324 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
1325
1326 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
1327 Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraCrsMat));
1328 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1329
1330 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1331 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
1332 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1333 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1334 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1335 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1336 return xpMat;
1337#else
1338 throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1339#endif
1340 }
1341#endif
1342
1343#ifdef HAVE_XPETRA_EPETRA
1344 if (isEpetra(op)) {
1345 Teuchos::RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*op);
1346 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1347 Teuchos::RCP<Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetra_op, true);
1348 Teuchos::RCP<Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetra_rowmat, true);
1349
1350 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpetraCrsMat =
1351 Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>(epetra_crsmat));
1352 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1353
1354 Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1355 Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xEpetraCrsMat, true);
1356 Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1357 Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1358 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1359 Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1360 return xpMat;
1361 }
1362#endif
1363 return Teuchos::null;
1364 }
1365
1366 static Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1367 toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1368#ifdef HAVE_XPETRA_TPETRA
1369 if (isTpetra(op)) {
1370 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1371 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
1372
1373 Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
1374 Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
1375 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
1376
1377 Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
1378 Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
1379 return xpOp;
1380 }
1381#endif
1382
1383#ifdef HAVE_XPETRA_EPETRA
1384 if (isEpetra(op)) {
1385 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
1386 }
1387#endif
1388 return Teuchos::null;
1389 }
1390
1391 static Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1392 toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
1393#ifdef HAVE_XPETRA_TPETRA
1394 if (isTpetra(op)) {
1395 typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1396 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
1397
1398 Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
1399 Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
1400 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
1401
1402 Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
1403 Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
1404 return xpOp;
1405 }
1406#endif
1407
1408#ifdef HAVE_XPETRA_EPETRA
1409 if (isEpetra(op)) {
1410 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
1411 }
1412#endif
1413 return Teuchos::null;
1414 }
1415
1416 static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1417 toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op) {
1418 using Teuchos::rcp_const_cast;
1419 using Teuchos::rcp_dynamic_cast;
1420 using ThyVSBase = Thyra::SpmdVectorSpaceBase<Scalar>;
1421 using thyTpV = Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
1422 using tV = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
1423
1424 RCP<const Thyra::VectorBase<Scalar>> diag = op->getDiag();
1425
1426 RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpDiag;
1427#ifdef HAVE_XPETRA_TPETRA
1428 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
1429 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
1430 if (!tDiag.is_null())
1431 xpDiag = Xpetra::toXpetra(tDiag);
1432 }
1433#endif
1434#ifdef HAVE_XPETRA_EPETRA
1435 if (xpDiag.is_null()) {
1436 RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(op->range())->getComm());
1437 RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(op->range()), comm);
1438 if (!map.is_null()) {
1439 RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
1440 RCP<Epetra_Vector> nceDiag = rcp_const_cast<Epetra_Vector>(eDiag);
1441 RCP<Xpetra::EpetraVectorT<int, Node>> xpEpDiag = rcp(new Xpetra::EpetraVectorT<int, Node>(nceDiag));
1442 xpDiag = rcp_dynamic_cast<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpEpDiag, true);
1443 }
1444 }
1445#endif
1446 TEUCHOS_ASSERT(!xpDiag.is_null());
1447 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
1448 return M;
1449 }
1450
1451 static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1452 toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op) {
1453 return toXpetra(Teuchos::rcp_const_cast<Thyra::DiagonalLinearOpBase<Scalar>>(op));
1454 }
1455
1456 static Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
1457 toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map) {
1458 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> thyraMap = Teuchos::null;
1459
1460 // check whether map is of type BlockedMap
1461 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
1462 if (bmap.is_null() == false) {
1463 Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>> vecSpaces(bmap->getNumMaps());
1464 for (size_t i = 0; i < bmap->getNumMaps(); i++) {
1465 // we need Thyra GIDs for all the submaps
1466 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> vs =
1467 Xpetra::ThyraUtils<Scalar, LO, GO, Node>::toThyra(bmap->getMap(i, true));
1468 vecSpaces[i] = vs;
1469 }
1470
1471 thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
1472 return thyraMap;
1473 }
1474
1475 // standard case
1476#ifdef HAVE_XPETRA_TPETRA
1477 if (map->lib() == Xpetra::UseTpetra) {
1478#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1479 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1480 Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>> tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(map);
1481 if (tpetraMap == Teuchos::null)
1482 throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
1483 RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tpMap = tpetraMap->getTpetra_Map();
1484 RCP<Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
1485 thyraMap = thyraTpetraMap;
1486#else
1487 throw Xpetra::Exceptions::RuntimeError("Problem DDD. Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1488#endif
1489 }
1490#endif
1491
1492#ifdef HAVE_XPETRA_EPETRA
1493 if (map->lib() == Xpetra::UseEpetra) {
1494 Teuchos::RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>> epetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(map);
1495 if (epetraMap == Teuchos::null)
1496 throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
1497 RCP<const Thyra::VectorSpaceBase<Scalar>> thyraEpetraMap = Thyra::create_VectorSpace(epetraMap->getEpetra_MapRCP());
1498 thyraMap = thyraEpetraMap;
1499 }
1500#endif
1501
1502 return thyraMap;
1503 }
1504
1505 static Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
1506 toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
1507 // create Thyra MultiVector
1508#ifdef HAVE_XPETRA_TPETRA
1509 if (vec->getMap()->lib() == Xpetra::UseTpetra) {
1510 auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
1511 RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpMV = Teuchos::rcp_dynamic_cast<const TpetraMultiVector>(vec)->getTpetra_MultiVector();
1512 auto thyDomMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(vec->getNumVectors(), vec->getMap()->getComm()));
1513 auto thyMV = rcp(new Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
1514 thyMV->initialize(thyTpMap, thyDomMap, tpMV);
1515 return thyMV;
1516 }
1517#endif
1518
1519#ifdef HAVE_XPETRA_EPETRA
1520 if (vec->getMap()->lib() == Xpetra::UseEpetra) {
1521 auto thyEpMap = Thyra::create_VectorSpace(Teuchos::rcp_dynamic_cast<const EpetraMapT<GlobalOrdinal, Node>>(vec->getMap())->getEpetra_MapRCP());
1522 auto epMV = Teuchos::rcp_dynamic_cast<const EpetraMultiVectorT<GlobalOrdinal, Node>>(vec)->getEpetra_MultiVector();
1523 auto thyMV = Thyra::create_MultiVector(epMV, thyEpMap);
1524 return thyMV;
1525 }
1526#endif
1527
1528 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "MultiVector cannot be converted to Thyra.");
1529 }
1530
1531 static Teuchos::RCP<const Thyra::VectorBase<Scalar>>
1532 toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
1533 // create Thyra Vector
1534#ifdef HAVE_XPETRA_TPETRA
1535 if (vec->getMap()->lib() == Xpetra::UseTpetra) {
1536 auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
1537 RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpVec = Teuchos::rcp_dynamic_cast<const TpetraVector>(vec)->getTpetra_Vector();
1538 auto thyVec = rcp(new Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
1539 thyVec->initialize(thyTpMap, tpVec);
1540 return thyVec;
1541 }
1542#endif
1543
1544#ifdef HAVE_XPETRA_EPETRA
1545 if (vec->getMap()->lib() == Xpetra::UseEpetra) {
1546 auto thyEpMap = Thyra::create_VectorSpace(Teuchos::rcp_dynamic_cast<const EpetraMapT<GlobalOrdinal, Node>>(vec->getMap())->getEpetra_MapRCP());
1547 auto epVec = rcp(Teuchos::rcp_dynamic_cast<const EpetraVectorT<GlobalOrdinal, Node>>(vec)->getEpetra_Vector(), false);
1548 auto thyVec = Thyra::create_Vector(epVec, thyEpMap);
1549 return thyVec;
1550 }
1551#endif
1552
1553 TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Vector cannot be converted to Thyra.");
1554 }
1555
1556 static void 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) {
1557 using Teuchos::as;
1558 using Teuchos::RCP;
1559 using Teuchos::rcp_dynamic_cast;
1560 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1561 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
1562 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1563 // typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
1564 // typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
1565 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1566
1567 // copy data from tY_inout to Y_inout
1568 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
1569 if (prodTarget != Teuchos::null) {
1570 RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
1571 if (bSourceVec.is_null() == true) {
1572 // SPECIAL CASE: target vector is product vector:
1573 // update Thyra product multi vector with data from a merged Xpetra multi vector
1574
1575 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
1576 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.");
1577
1578 for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1579 // access Xpetra data
1580 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
1581
1582 // access Thyra data
1583 Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
1584 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
1585 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
1586 const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
1587 const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim());
1588 RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
1589 Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
1590
1591 // loop over all vectors in multivector
1592 for (size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
1593 Teuchos::ArrayRCP<const Scalar> xpData = xpSubBlock->getData(j); // access const data from Xpetra object
1594
1595 // loop over all local rows
1596 for (LocalOrdinal i = 0; i < localSubDim; ++i) {
1597 (*thyData)(i, j) = xpData[i];
1598 }
1599 }
1600 }
1601 } else {
1602 // source vector is a blocked multivector
1603 // TODO test me
1604 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.");
1605
1606 for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1607 // access Thyra data
1608 RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
1609
1610 Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
1611
1612 // access Thyra data
1613 Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
1614 Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
1615 }
1616 }
1617 } else {
1618 // STANDARD case:
1619 // update Thyra::MultiVector with data from an Xpetra::MultiVector
1620
1621 // access Thyra data
1622 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
1623 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
1624 const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
1625 const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim());
1626 RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
1627 Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
1628
1629 // loop over all vectors in multivector
1630 for (size_t j = 0; j < source->getNumVectors(); ++j) {
1631 Teuchos::ArrayRCP<const Scalar> xpData = source->getData(j); // access const data from Xpetra object
1632 // loop over all local rows
1633 for (LocalOrdinal i = 0; i < localSubDim; ++i) {
1634 (*thyData)(i, j) = xpData[i];
1635 }
1636 }
1637 }
1638 }
1639
1640 static Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
1641 toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
1642 // create a Thyra operator from Xpetra::CrsMatrix
1643 Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
1644
1645#ifdef HAVE_XPETRA_TPETRA
1646 Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
1647 if (tpetraMat != Teuchos::null) {
1648#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1649 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1650
1651 Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
1652 Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
1653 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
1654
1655 Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
1656 Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
1657
1658 thyraOp = Thyra::createConstLinearOp(tpOperator);
1659 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
1660#else
1661 throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1662#endif
1663 }
1664#endif
1665
1666#ifdef HAVE_XPETRA_EPETRA
1667 Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
1668 if (epetraMat != Teuchos::null) {
1669 Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat, true);
1670 Teuchos::RCP<const Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrix();
1671 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
1672
1673 Teuchos::RCP<const Thyra::EpetraLinearOp> thyraEpOp = Thyra::epetraLinearOp(epCrsMat, "op");
1674 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
1675 thyraOp = thyraEpOp;
1676 }
1677#endif
1678 return thyraOp;
1679 }
1680
1681 static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
1682 toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
1683 // create a Thyra operator from Xpetra::CrsMatrix
1684 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
1685
1686#ifdef HAVE_XPETRA_TPETRA
1687 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
1688 if (tpetraMat != Teuchos::null) {
1689#if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1690 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1691
1692 Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
1693 Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
1694 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
1695
1696 Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
1697 Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
1698
1699 thyraOp = Thyra::createLinearOp(tpOperator);
1700 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
1701#else
1702 throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1703#endif
1704 }
1705#endif
1706
1707#ifdef HAVE_XPETRA_EPETRA
1708 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
1709 if (epetraMat != Teuchos::null) {
1710 Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat, true);
1711 Teuchos::RCP<Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrixNonConst();
1712 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
1713
1714 Teuchos::RCP<Thyra::EpetraLinearOp> thyraEpOp = Thyra::nonconstEpetraLinearOp(epCrsMat, "op");
1715 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
1716 thyraOp = thyraEpOp;
1717 }
1718#endif
1719 return thyraOp;
1720 }
1721
1722 static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
1723 toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, long long, EpetraNode>>& mat);
1724
1725}; // specialization on SC=double, LO=GO=int and NO=EpetraNode
1726#endif // XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
1727
1728#endif // HAVE_XPETRA_EPETRA
1729
1730} // end namespace Xpetra
1731
1732#define XPETRA_THYRAUTILS_SHORT
1733#endif // HAVE_XPETRA_THYRA
1734
1735#endif // XPETRA_THYRAUTILS_HPP
Xpetra_ConfigDefs.hpp
EpetraNode
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
Definition Xpetra_ETI_3arg.hpp:31
Xpetra_EpetraCrsMatrix.hpp
Xpetra_EpetraMap.hpp
Xpetra_MapUtils.hpp
Xpetra_TpetraOperator.hpp
Xpetra_UseShortNames.hpp
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::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition Xpetra_MatrixFactory_def.hpp:22
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::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