Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
76 namespace Xpetra {
77 
78 template <class Scalar,
79  class LocalOrdinal = int,
80  class GlobalOrdinal = LocalOrdinal,
81  class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
82 class 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
156 template <>
157 class 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
942 template <>
943 class 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_TpetraOperator.hpp
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
Definition: Xpetra_MultiVectorFactory_def.hpp:24
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrixNonConst
RCP< Epetra_CrsMatrix > getEpetra_CrsMatrixNonConst() const
Definition: Xpetra_EpetraCrsMatrix.hpp:208
Xpetra_UseShortNames.hpp
Xpetra_MapUtils.hpp
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:65
Xpetra::MapUtils::transformThyra2XpetraGIDs
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > transformThyra2XpetraGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlInput, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlReferenceInput)
replace set of global ids by new global ids
Definition: Xpetra_MapUtils.hpp:191
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrix
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
Definition: Xpetra_EpetraCrsMatrix.hpp:207
Xpetra_EpetraCrsMatrix.hpp
Xpetra_EpetraMap.hpp
EpetraNode
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
Definition: Xpetra_ETI_3arg.hpp:31
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory_def.hpp:22
Xpetra::TpetraCrsMatrix::getTpetra_CrsMatrix
RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrix() const
Get the underlying Tpetra matrix.
Definition: Xpetra_TpetraCrsMatrix_def.hpp:611
Xpetra::EpetraMapT::getEpetra_MapRCP
const RCP< const Epetra_Map > & getEpetra_MapRCP() const
Get the underlying Epetra map.
Definition: Xpetra_EpetraMap.hpp:212
Xpetra_ConfigDefs.hpp
Xpetra::TpetraCrsMatrix::getTpetra_CrsMatrixNonConst
RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrixNonConst() const
Get the underlying Tpetra matrix.
Definition: Xpetra_TpetraCrsMatrix_def.hpp:615
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_decl.hpp:46
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
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_decl.hpp:60
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