Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Xpetra_ThyraUtils.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef XPETRA_THYRAUTILS_HPP
48 #define XPETRA_THYRAUTILS_HPP
49 
50 #include "Xpetra_ConfigDefs.hpp"
51 #ifdef HAVE_XPETRA_THYRA
52 
53 #include <typeinfo>
54 
55 #ifdef HAVE_XPETRA_TPETRA
56 #include "Tpetra_ConfigDefs.hpp"
57 #endif
58 
59 #ifdef HAVE_XPETRA_EPETRA
60 #include "Epetra_config.h"
61 #include "Epetra_CombineMode.h"
62 #endif
63 
64 #include "Xpetra_Map.hpp"
65 #include "Xpetra_BlockedMap.hpp"
66 #include "Xpetra_BlockedMultiVector.hpp"
67 #include "Xpetra_BlockedCrsMatrix.hpp"
68 #include "Xpetra_MapUtils.hpp"
69 #include "Xpetra_StridedMap.hpp"
70 #include "Xpetra_StridedMapFactory.hpp"
71 #include "Xpetra_MapExtractor.hpp"
72 #include "Xpetra_Matrix.hpp"
73 #include "Xpetra_MatrixFactory.hpp"
74 #include "Xpetra_CrsMatrixWrap.hpp"
75 #include "Xpetra_MultiVectorFactory.hpp"
76 
77 #include <Thyra_VectorSpaceBase.hpp>
78 #include <Thyra_SpmdVectorSpaceBase.hpp>
79 #include <Thyra_ProductVectorSpaceBase.hpp>
80 #include <Thyra_ProductMultiVectorBase.hpp>
81 #include <Thyra_VectorSpaceBase.hpp>
82 #include <Thyra_DefaultProductVectorSpace.hpp>
83 #include <Thyra_DefaultBlockedLinearOp.hpp>
84 #include <Thyra_LinearOpBase.hpp>
85 #include "Thyra_DiagonalLinearOpBase.hpp"
86 #include <Thyra_DetachedMultiVectorView.hpp>
87 #include <Thyra_MultiVectorStdOps.hpp>
88 
89 #ifdef HAVE_XPETRA_TPETRA
90 #include <Thyra_TpetraThyraWrappers.hpp>
91 #include <Thyra_TpetraVector.hpp>
92 #include <Thyra_TpetraMultiVector.hpp>
93 #include <Thyra_TpetraVectorSpace.hpp>
94 #include <Tpetra_Map.hpp>
95 #include <Tpetra_Vector.hpp>
96 #include <Tpetra_CrsMatrix.hpp>
97 #include <Xpetra_TpetraMap.hpp>
98 #include <Xpetra_TpetraOperator.hpp>
99 #include <Xpetra_TpetraCrsMatrix.hpp>
100 #endif
101 #ifdef HAVE_XPETRA_EPETRA
102 #include <Thyra_EpetraLinearOp.hpp>
103 #include <Thyra_EpetraThyraWrappers.hpp>
104 #include <Thyra_SpmdVectorBase.hpp>
105 #include <Thyra_get_Epetra_Operator.hpp>
106 #include <Epetra_Map.h>
107 #include <Epetra_Vector.h>
108 #include <Epetra_CrsMatrix.h>
109 #include <Xpetra_EpetraMap.hpp>
110 #include <Xpetra_EpetraCrsMatrix.hpp>
111 #endif
112 
113 namespace Xpetra {
114 
115 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
116 class BlockedCrsMatrix;
117 
118 template <class Scalar,
119  class LocalOrdinal = int,
120  class GlobalOrdinal = LocalOrdinal,
121  class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
122 class ThyraUtils {
123  private:
124 #undef XPETRA_THYRAUTILS_SHORT
125 #include "Xpetra_UseShortNames.hpp"
126 
127  public:
128  static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>
129  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) {
130  Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map = toXpetra(vectorSpace);
131 
132  if (stridedBlockId == -1) {
133  TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo.size() != 0);
134  } else {
135  TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo[stridedBlockId] != 0);
136  }
137 
138  Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>> ret = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(map, stridingInfo, stridedBlockId, offset);
139  return ret;
140  }
141 
142  static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
143  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
144  using Teuchos::as;
145  using Teuchos::RCP;
146  using Teuchos::rcp_dynamic_cast;
147  using ThyVecSpaceBase = Thyra::VectorSpaceBase<Scalar>;
148  using ThyProdVecSpaceBase = Thyra::ProductVectorSpaceBase<Scalar>;
149  using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
150 
151  RCP<Map> resultMap = Teuchos::null;
152  RCP<const ThyProdVecSpaceBase> prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase>(vectorSpace);
153  if (prodVectorSpace != Teuchos::null) {
154  // SPECIAL CASE: product Vector space
155  // collect all submaps to store them in a hierarchical BlockedMap object
156  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks() == 0, std::logic_error, "Found a product vector space with zero blocks.");
157  std::vector<RCP<const Map>> mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
158  std::vector<RCP<const Map>> mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
159  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
160  RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
161  // can be of type Map or BlockedMap (containing Thyra GIDs)
162  mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
163  }
164 
165  // get offsets for submap GIDs
166  // we need that for the full map (Xpetra GIDs)
167  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(), 0);
168  for (int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
169  gidOffsets[i] = mapsThyra[i - 1]->getMaxAllGlobalIndex() + gidOffsets[i - 1] + 1;
170  }
171 
172  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
173  RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
174  // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
175  mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
176  }
177 
178  resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapsXpetra, mapsThyra));
179  } else {
180 #ifdef HAVE_XPETRA_TPETRA
181  // STANDARD CASE: no product map
182  // check whether we have a Tpetra based Thyra operator
183  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vectorSpace);
184  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(tpetra_vsc) == true, Xpetra::Exceptions::RuntimeError, "toXpetra failed to convert provided vector space to Thyra::TpetraVectorSpace. This is the general implementation for Tpetra only.");
185 
186  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> TpMap;
187  typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpVector;
188  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
189  typedef Thyra::VectorBase<Scalar> ThyVecBase;
190  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
191  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
192  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
193  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
194  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
195  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
196  resultMap = Xpetra::toXpetraNonConst(rgTpetraMap);
197 #else
198  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Cannot transform Thyra::VectorSpace to Xpetra::Map. This is the general implementation for Tpetra only, but Tpetra is disabled.");
199 #endif
200  } // end standard case (no product map)
201  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
202  return resultMap;
203  }
204 
205  // const version
206  static Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
207  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
208  using Teuchos::as;
209  using Teuchos::RCP;
210  using Teuchos::rcp_dynamic_cast;
211 
212  using ThyMultVecBase = Thyra::MultiVectorBase<Scalar>;
213  using ThyProdMultVecBase = Thyra::ProductMultiVectorBase<Scalar>;
214  using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
215 
216  // return value
217  RCP<MultiVector> xpMultVec = Teuchos::null;
218 
219  // check whether v is a product multi vector
220  Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase>(v);
221  if (thyProdVec != Teuchos::null) {
222  // SPECIAL CASE: create a nested BlockedMultiVector
223  // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
224  RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
225  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
226 
227  // create new Xpetra::BlockedMultiVector
228  xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
229 
230  RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec, true);
231 
232  // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
233  for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b) {
234  RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
235  // xpBlockMV can be of type MultiVector or BlockedMultiVector
236  RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); // recursive call
237  xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
238  }
239  } else {
240  // STANDARD CASE: no product vector
241 #ifdef HAVE_XPETRA_TPETRA
242  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> ConverterT;
243  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpMultVec;
244  typedef Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpTpMultVec;
245  typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
246  typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
247 
248  RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
249  RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
250  TEUCHOS_TEST_FOR_EXCEPTION(thyraTpetraMultiVector == Teuchos::null, Xpetra::Exceptions::RuntimeError, "toXpetra failed to convert provided multi vector to Thyra::TpetraMultiVector. This is the general implementation for Tpetra only.");
251  const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
252  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
253  RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
254  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
255  xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
256 #else
257  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Cannot transform Thyra::MultiVector to Xpetra::MultiVector. This is the general implementation for Tpetra only, but Teptra is disabled.");
258 #endif
259  } // end standard case
260  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
261  return xpMultVec;
262  }
263 
264  // non-const version
265  static Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
266  toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
267  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> cv =
268  Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar>>(v);
269  Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> r =
270  toXpetra(cv, comm);
271  return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(r);
272  }
273 
274  static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
275  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(op));
276 
277  // check whether we have a Tpetra based Thyra operator
278  bool bIsTpetra = false;
279 #ifdef HAVE_XPETRA_TPETRA
280  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(op);
281  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
282 
283  // for debugging purposes: find out why dynamic cast failed
284  if (!bIsTpetra &&
285 #ifdef HAVE_XPETRA_EPETRA
286  Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
287 #endif
288  Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar>>(op) == Teuchos::null) {
289  // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
290  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
291  if (Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
292  std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
293  std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
294  std::cout << " properly set!" << std::endl;
295  std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
296  }
297  }
298 #endif
299 
300 #if 0
301  // Check whether it is a blocked operator.
302  // If yes, grab the (0,0) block and check the underlying linear algebra
303  // Note: we expect that the (0,0) block exists!
304  if(bIsTpetra == false) {
305  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
306  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
307  if(ThyBlockedOp != Teuchos::null) {
308  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
309  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
310  ThyBlockedOp->getBlock(0,0);
311  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
312  bIsTpetra = isTpetra(b00);
313  }
314  }
315 #endif
316 
317  return bIsTpetra;
318  }
319 
320  static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
321  return false;
322  }
323 
324  static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
325  // Check whether it is a blocked operator.
326  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar>> ThyBlockedOp =
327  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar>>(op);
328  if (ThyBlockedOp != Teuchos::null) {
329  return true;
330  }
331  return false;
332  }
333 
334  static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
335  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
336 #ifdef HAVE_XPETRA_TPETRA
337  if (isTpetra(op)) {
338  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
339  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
340  // we should also add support for the const versions!
341  // getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
342  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
343  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
344  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
345  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat);
346  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
347 
348  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
349  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraNcnstCrsMat));
350  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
351 
352  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
353  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
354  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
355  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
356  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
357  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
358  return xpMat;
359  }
360 #endif
361 
362 #ifdef HAVE_XPETRA_EPETRA
363  if (isEpetra(op)) {
364  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
365  }
366 #endif
367  return Teuchos::null;
368  }
369 
370  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
371  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
372  using Teuchos::RCP;
373  using Teuchos::rcp;
374  using Teuchos::rcp_const_cast;
375  using Teuchos::rcp_dynamic_cast;
376 
377 #ifdef HAVE_XPETRA_TPETRA
378  if (isTpetra(op)) {
379  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
380  typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraOperator_t;
381  RCP<const TpetraOperator_t> TpetraOp = TOE::getConstTpetraOperator(op);
382  TEUCHOS_TEST_FOR_EXCEPT(is_null(TpetraOp));
383  RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat =
384  rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
385  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat =
386  rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
387 
388  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
389  rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(
390  rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat)));
391  TEUCHOS_TEST_FOR_EXCEPT(is_null(xTpetraCrsMat));
392  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
393  rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
394  RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
395  rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
396  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
397  rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
398  return xpMat;
399  }
400 #endif
401 
402 #ifdef HAVE_XPETRA_EPETRA
403  if (isEpetra(op)) {
404  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
405  }
406 #endif
407  return Teuchos::null;
408  }
409 
410  static Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
411  toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
412 #ifdef HAVE_XPETRA_TPETRA
413  if (isTpetra(op)) {
414  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
415  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
416 
417  Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
418  Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
419  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
420 
421  Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
422  Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
423  return xpOp;
424  }
425 #endif
426 
427 #ifdef HAVE_XPETRA_EPETRA
428  if (isEpetra(op)) {
429  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
430  }
431 #endif
432  return Teuchos::null;
433  }
434 
435  static Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
436  toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
437 #ifdef HAVE_XPETRA_TPETRA
438  if (isTpetra(op)) {
439  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
440  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
441 
442  Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
443  Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
444  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
445 
446  Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
447  Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
448  return xpOp;
449  }
450 #endif
451 
452 #ifdef HAVE_XPETRA_EPETRA
453  if (isEpetra(op)) {
454  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
455  }
456 #endif
457  return Teuchos::null;
458  }
459 
460  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
461  toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op) {
462  using Teuchos::rcp_const_cast;
463  using Teuchos::rcp_dynamic_cast;
464 
465  RCP<const Thyra::VectorBase<Scalar>> diag = op->getDiag();
466 
467  RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpDiag;
468 #ifdef HAVE_XPETRA_TPETRA
469  using thyTpV = Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
470  using tV = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
471  if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
472  RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
473  if (!tDiag.is_null())
474  xpDiag = Xpetra::toXpetra(tDiag);
475  }
476 #endif
477  TEUCHOS_ASSERT(!xpDiag.is_null());
478  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
479  return M;
480  }
481 
482  static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
483  toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op) {
484  return toXpetra(Teuchos::rcp_const_cast<Thyra::DiagonalLinearOpBase<Scalar>>(op));
485  }
486 
487  static Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
488  toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map) {
489  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> thyraMap = Teuchos::null;
490 
491  // check whether map is of type BlockedMap
492  RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
493  if (bmap.is_null() == false) {
494  Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>> vecSpaces(bmap->getNumMaps());
495  for (size_t i = 0; i < bmap->getNumMaps(); i++) {
496  // we need Thyra GIDs for all the submaps
497  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> vs =
498  Xpetra::ThyraUtils<Scalar, LO, GO, Node>::toThyra(bmap->getMap(i, true));
499  vecSpaces[i] = vs;
500  }
501 
502  thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
503  return thyraMap;
504  }
505 
506  // standard case
507 #ifdef HAVE_XPETRA_TPETRA
508  if (map->lib() == Xpetra::UseTpetra) {
509  Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>> tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(map);
510  if (tpetraMap == Teuchos::null)
511  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
512  RCP<Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetraMap->getTpetra_Map());
513  thyraMap = thyraTpetraMap;
514  }
515 #endif
516 
517 #ifdef HAVE_XPETRA_EPETRA
518  if (map->lib() == Xpetra::UseEpetra) {
519  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
520  }
521 #endif
522 
523  return thyraMap;
524  }
525 
526  static Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
527  toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
528  // create Thyra MultiVector
529 #ifdef HAVE_XPETRA_TPETRA
530  if (vec->getMap()->lib() == Xpetra::UseTpetra) {
531  auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
532  RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpMV = Teuchos::rcp_dynamic_cast<const TpetraMultiVector>(vec)->getTpetra_MultiVector();
533  auto thyDomMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(vec->getNumVectors(), vec->getMap()->getComm()));
534  auto thyMV = rcp(new Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
535  thyMV->initialize(thyTpMap, thyDomMap, tpMV);
536  return thyMV;
537  }
538 #endif
539 
540 #ifdef HAVE_XPETRA_EPETRA
541  if (vec->getMap()->lib() == Xpetra::UseEpetra) {
542  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
543  }
544 #endif
545 
546  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "MultiVector cannot be converted to Thyra.");
547  }
548 
549  static Teuchos::RCP<const Thyra::VectorBase<Scalar>>
550  toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
551  // create Thyra Vector
552 #ifdef HAVE_XPETRA_TPETRA
553  if (vec->getMap()->lib() == Xpetra::UseTpetra) {
554  auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
555  RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpVec = Teuchos::rcp_dynamic_cast<const TpetraVector>(vec)->getTpetra_Vector();
556  auto thyVec = rcp(new Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
557  thyVec->initialize(thyTpMap, tpVec);
558  return thyVec;
559  }
560 #endif
561 
562 #ifdef HAVE_XPETRA_EPETRA
563  if (vec->getMap()->lib() == Xpetra::UseEpetra) {
564  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
565  }
566 #endif
567 
568  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Vector cannot be converted to Thyra.");
569  }
570 
571  // update Thyra multi vector with data from Xpetra multi vector
572  // In case of a Thyra::ProductMultiVector the Xpetra::MultiVector is splitted into its subparts using a provided MapExtractor
573  static void
574  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) {
575  using Teuchos::as;
576  using Teuchos::RCP;
577  using Teuchos::rcp_dynamic_cast;
578  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
579  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
580  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
581  // typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
582  // typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
583  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
584 
585  // copy data from tY_inout to Y_inout
586  RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
587  if (prodTarget != Teuchos::null) {
588  RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
589  if (bSourceVec.is_null() == true) {
590  // SPECIAL CASE: target vector is product vector:
591  // update Thyra product multi vector with data from a merged Xpetra multi vector
592 
593  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
594  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.");
595 
596  for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
597  // access Xpetra data
598  RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
599 
600  // access Thyra data
601  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
602  RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
603  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
604  const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
605  const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim());
606  RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
607  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
608 
609  // loop over all vectors in multivector
610  for (size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
611  Teuchos::ArrayRCP<const Scalar> xpData = xpSubBlock->getData(j); // access const data from Xpetra object
612 
613  // loop over all local rows
614  for (LocalOrdinal i = 0; i < localSubDim; ++i) {
615  (*thyData)(i, j) = xpData[i];
616  }
617  }
618  }
619  } else {
620  // source vector is a blocked multivector
621  // TODO test me
622  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.");
623 
624  for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
625  // access Thyra data
626  RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
627 
628  Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
629 
630  // access Thyra data
631  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
632  Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
633  }
634  }
635  } else {
636  // STANDARD case:
637  // update Thyra::MultiVector with data from an Xpetra::MultiVector
638 
639  // access Thyra data
640  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
641  TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
642  const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
643  const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim());
644  RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
645  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
646 
647  // loop over all vectors in multivector
648  for (size_t j = 0; j < source->getNumVectors(); ++j) {
649  Teuchos::ArrayRCP<const Scalar> xpData = source->getData(j); // access const data from Xpetra object
650  // loop over all local rows
651  for (LocalOrdinal i = 0; i < localSubDim; ++i) {
652  (*thyData)(i, j) = xpData[i];
653  }
654  }
655  }
656  }
657 
658  static Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
659  toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
660  // create a Thyra operator from Xpetra::CrsMatrix
661  Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
662 
663  // bool bIsTpetra = false;
664 
665 #ifdef HAVE_XPETRA_TPETRA
666  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
667  if (tpetraMat != Teuchos::null) {
668  // bIsTpetra = true;
669  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
670  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
671  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
672 
673  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
674  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
675 
676  thyraOp = Thyra::createConstLinearOp(tpOperator);
677  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
678  } else {
679 #ifdef HAVE_XPETRA_EPETRA
680  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to Tpetra::CrsMatrix failed. Assume matrix should be Epetra then. Epetra needs SC=double, LO=int, and GO=int or GO=long long");
681 #else
682  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to Tpetra::CrsMatrix failed. Assume matrix should be Epetra then. No Epetra available");
683 #endif
684  }
685  return thyraOp;
686 #else
687  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
688  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
689 #endif
690  }
691 
692  static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
693  toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
694  // create a Thyra operator from Xpetra::CrsMatrix
695  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
696 
697  // bool bIsTpetra = false;
698 
699 #ifdef HAVE_XPETRA_TPETRA
700  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
701  if (tpetraMat != Teuchos::null) {
702  // bIsTpetra = true;
703  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
704  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
705  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
706 
707  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
708  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
709 
710  thyraOp = Thyra::createLinearOp(tpOperator);
711  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
712  } else {
713  // cast to TpetraCrsMatrix failed
714 #ifdef HAVE_XPETRA_EPETRA
715  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Assuming matrix supposed to be Epetra. Epetra needs SC=double, LO=int, and GO=int or GO=long long");
716 #else
717  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Guess, matrix should be Epetra then, but no Epetra available.");
718 #endif
719  }
720  return thyraOp;
721 #else
722  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
723  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
724 #endif
725  }
726 
727  static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
728  toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
729  int nRows = mat->Rows();
730  int nCols = mat->Cols();
731 
732  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Ablock = mat->getInnermostCrsMatrix();
733  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock);
734  TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
735 
736 #ifdef HAVE_XPETRA_TPETRA
737  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
738  if (tpetraMat != Teuchos::null) {
739  // create new Thyra blocked operator
740  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar>> blockMat =
741  Thyra::defaultBlockedLinearOp<Scalar>();
742 
743  blockMat->beginBlockFill(nRows, nCols);
744 
745  for (int r = 0; r < nRows; ++r) {
746  for (int c = 0; c < nCols; ++c) {
747  Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r, c);
748 
749  if (xpmat == Teuchos::null) continue; // shortcut for empty blocks
750 
751  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thBlock = Teuchos::null;
752 
753  // check whether the subblock is again a blocked operator
754  Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpblock =
755  Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpmat);
756  if (xpblock != Teuchos::null) {
757  if (xpblock->Rows() == 1 && xpblock->Cols() == 1) {
758  // If it is a single block operator, unwrap it
759  Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
760  TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
761  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
762  } else {
763  // recursive call for general blocked operators
764  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpblock);
765  }
766  } else {
767  // check whether it is a CRSMatrix object
768  Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
769  TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
770  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
771  }
772 
773  blockMat->setBlock(r, c, thBlock);
774  }
775  }
776 
777  blockMat->endBlockFill();
778 
779  return blockMat;
780  } else {
781  // tpetraMat == Teuchos::null
782 #ifdef HAVE_XPETRA_EPETRA
783  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Assuming matrix supposed to be Epetra. Epetra needs SC=double, LO=int, and GO=int or GO=long long");
784 #else
785  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Guess, matrix should be Epetra then, but no Epetra available.");
786 #endif
787  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
788  }
789 #endif // endif HAVE_XPETRA_TPETRA
790 
791 // 4-Aug-2017 JJH Added 2nd condition to avoid "warning: dynamic initialization in unreachable code"
792 // If HAVE_XPETRA_TPETRA is defined, then this method will always return or throw in the if-then-else above.
793 #if defined(HAVE_XPETRA_EPETRA) && !defined(HAVE_XPETRA_TPETRA)
794  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
795  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
796 #endif // endif HAVE_XPETRA_EPETRA
797  }
798 
799 }; // end Utils class
800 
801 // full specialization for Epetra support
802 // Note, that Thyra only has support for Epetra (not for Epetra64)
803 #ifdef HAVE_XPETRA_EPETRA
804 
805 #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
806 template <>
807 class ThyraUtils<double, int, int, EpetraNode> {
808  public:
809  typedef double Scalar;
810  typedef int LocalOrdinal;
811  typedef int GlobalOrdinal;
812  typedef EpetraNode Node;
813 
814  private:
815 #undef XPETRA_THYRAUTILS_SHORT
816 #include "Xpetra_UseShortNames.hpp"
817 
818  public:
819  static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>
820  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) {
821  Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map = ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(vectorSpace, comm);
822 
823  if (stridedBlockId == -1) {
824  TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo.size() != 0);
825  } else {
826  TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo[stridedBlockId] != 0);
827  }
828 
829  Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>> ret = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(map, stridingInfo, stridedBlockId, offset);
830  return ret;
831  }
832 
833  static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
834  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
835  using Teuchos::as;
836  using Teuchos::RCP;
837  using Teuchos::rcp_dynamic_cast;
838  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
839  typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
840  typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyUtils;
841 
842  RCP<Map> resultMap = Teuchos::null;
843 
844  RCP<const ThyProdVecSpaceBase> prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase>(vectorSpace);
845  if (prodVectorSpace != Teuchos::null) {
846  // SPECIAL CASE: product Vector space
847  // collect all submaps to store them in a hierarchical BlockedMap object
848  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks() == 0, std::logic_error, "Found a product vector space with zero blocks.");
849  std::vector<RCP<const Map>> mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
850  std::vector<RCP<const Map>> mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
851  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
852  RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
853  // can be of type Map or BlockedMap (containing Thyra GIDs)
854  mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
855  }
856 
857  // get offsets for submap GIDs
858  // we need that for the full map (Xpetra GIDs)
859  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(), 0);
860  for (int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
861  gidOffsets[i] = mapsThyra[i - 1]->getMaxAllGlobalIndex() + gidOffsets[i - 1] + 1;
862  }
863 
864  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
865  RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
866  // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
867  mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
868  }
869 
870  resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapsXpetra, mapsThyra));
871  } else {
872  // STANDARD CASE: no product map
873  // Epetra/Tpetra specific code to access the underlying map data
874 
875  // check whether we have a Tpetra based Thyra operator
876  bool bIsTpetra = false;
877 #ifdef HAVE_XPETRA_TPETRA
878 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
879  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
880  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vectorSpace);
881  bIsTpetra = Teuchos::is_null(tpetra_vsc) ? false : true;
882 #endif
883 #endif
884 
885  // check whether we have an Epetra based Thyra operator
886  bool bIsEpetra = !bIsTpetra; // note: this is a little bit fragile!
887 
888 #ifdef HAVE_XPETRA_TPETRA
889  if (bIsTpetra) {
890 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
891  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
892  typedef Thyra::VectorBase<Scalar> ThyVecBase;
893  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> TpMap;
894  typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpVector;
895  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
896  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
897  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
898  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
899  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
900  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
901  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
902 
903  resultMap = Xpetra::toXpetraNonConst(rgTpetraMap);
904 #else
905  throw Xpetra::Exceptions::RuntimeError("Problem AAA. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
906 #endif
907  }
908 #endif
909 
910 #ifdef HAVE_XPETRA_EPETRA
911  if (bIsEpetra) {
912  // RCP<const Epetra_Map> epMap = Teuchos::null;
913  RCP<const Epetra_Map>
914  epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map>>(vectorSpace, "epetra_map");
915  if (!Teuchos::is_null(epetra_map)) {
916  resultMap = Teuchos::rcp(new Xpetra::EpetraMapT<GlobalOrdinal, Node>(epetra_map));
917  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
918  } else {
919  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "No Epetra_Map data found in Thyra::VectorSpace.");
920  }
921  }
922 #endif
923  } // end standard case (no product map)
924  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
925  return resultMap;
926  }
927 
928  // const version
929  static Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
930  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
931  using Teuchos::as;
932  using Teuchos::RCP;
933  using Teuchos::rcp_dynamic_cast;
934 
935  using ThyProdMultVecBase = Thyra::ProductMultiVectorBase<Scalar>;
936  using ThyMultVecBase = Thyra::MultiVectorBase<Scalar>;
937  using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
938 
939  // return value
940  RCP<MultiVector> xpMultVec = Teuchos::null;
941 
942  // check whether v is a product multi vector
943  Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase>(v);
944  if (thyProdVec != Teuchos::null) {
945  // SPECIAL CASE: create a nested BlockedMultiVector
946  // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
947  RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
948  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
949 
950  // create new Xpetra::BlockedMultiVector
951  xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
952 
953  RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec, true);
954 
955  // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
956  for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b) {
957  RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
958  // xpBlockMV can be of type MultiVector or BlockedMultiVector
959  RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); // recursive call
960  xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
961  }
962 
963  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
964  return xpMultVec;
965  } else {
966  // STANDARD CASE: no product vector
967  // Epetra/Tpetra specific code to access the underlying map data
968  bool bIsTpetra = false;
969 #ifdef HAVE_XPETRA_TPETRA
970 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
971  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
972 
973  // typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
974  // typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
975  typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
976  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> ConverterT;
977  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpMultVec;
978  typedef Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpTpMultVec;
979  typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
980 
981  RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
982  RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
983  if (thyraTpetraMultiVector != Teuchos::null) {
984  bIsTpetra = true;
985  const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
986  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
987  RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
988  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
989  xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
990  }
991 #endif
992 #endif
993 
994 #ifdef HAVE_XPETRA_EPETRA
995  if (bIsTpetra == false) {
996  // no product vector
997  Teuchos::RCP<Map> map = ThyUtils::toXpetra(v->range(), comm);
998  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(map));
999  RCP<Xpetra::EpetraMapT<GlobalOrdinal, Node>> xeMap = rcp_dynamic_cast<Xpetra::EpetraMapT<GlobalOrdinal, Node>>(map, true);
1000  RCP<const Epetra_Map> eMap = xeMap->getEpetra_MapRCP();
1001  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(eMap));
1002  Teuchos::RCP<const Epetra_MultiVector> epMultVec = Thyra::get_Epetra_MultiVector(*eMap, v);
1003  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epMultVec));
1004  RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<Epetra_MultiVector>(epMultVec);
1005  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epNonConstMultVec));
1006  xpMultVec = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>(epNonConstMultVec));
1007  }
1008 #endif
1009  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1010  return xpMultVec;
1011  } // end standard case
1012  }
1013 
1014  // non-const version
1015  static Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1016  toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
1017  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> cv =
1018  Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar>>(v);
1019  Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> r =
1020  toXpetra(cv, comm);
1021  return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(r);
1022  }
1023 
1024  static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1025  // check whether we have a Tpetra based Thyra operator
1026  bool bIsTpetra = false;
1027 #ifdef HAVE_XPETRA_TPETRA
1028 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1029  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1030 
1031  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(op);
1032  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
1033 
1034  // for debugging purposes: find out why dynamic cast failed
1035  if (!bIsTpetra &&
1036 #ifdef HAVE_XPETRA_EPETRA
1037  Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
1038 #endif
1039  Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar>>(op) == Teuchos::null) {
1040  // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
1041  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
1042  if (Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
1043  std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
1044  std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
1045  std::cout << " properly set!" << std::endl;
1046  std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
1047  }
1048  }
1049 #endif
1050 #endif
1051 
1052 #if 0
1053  // Check whether it is a blocked operator.
1054  // If yes, grab the (0,0) block and check the underlying linear algebra
1055  // Note: we expect that the (0,0) block exists!
1056  if(bIsTpetra == false) {
1057  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1058  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
1059  if(ThyBlockedOp != Teuchos::null) {
1060  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1061  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1062  ThyBlockedOp->getBlock(0,0);
1063  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1064  bIsTpetra = isTpetra(b00);
1065  }
1066  }
1067 #endif
1068 
1069  return bIsTpetra;
1070  }
1071 
1072  static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1073  // check whether we have an Epetra based Thyra operator
1074  bool bIsEpetra = false;
1075 
1076 #ifdef HAVE_XPETRA_EPETRA
1077  Teuchos::RCP<const Thyra::EpetraLinearOp> epetra_op = Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op, false);
1078  bIsEpetra = Teuchos::is_null(epetra_op) ? false : true;
1079 #endif
1080 
1081 #if 0
1082  // Check whether it is a blocked operator.
1083  // If yes, grab the (0,0) block and check the underlying linear algebra
1084  // Note: we expect that the (0,0) block exists!
1085  if(bIsEpetra == false) {
1086  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1087  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op,false);
1088  if(ThyBlockedOp != Teuchos::null) {
1089  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1090  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1091  ThyBlockedOp->getBlock(0,0);
1092  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1093  bIsEpetra = isEpetra(b00);
1094  }
1095  }
1096 #endif
1097 
1098  return bIsEpetra;
1099  }
1100 
1101  static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1102  // Check whether it is a blocked operator.
1103  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar>> ThyBlockedOp =
1104  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar>>(op);
1105  if (ThyBlockedOp != Teuchos::null) {
1106  return true;
1107  }
1108  return false;
1109  }
1110 
1111  static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1112  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1113 #ifdef HAVE_XPETRA_TPETRA
1114  if (isTpetra(op)) {
1115 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1116  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1117 
1118  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1119  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
1120  // we should also add support for the const versions!
1121  // getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
1122  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1123  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp);
1124  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraRowMat));
1125  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
1126  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat);
1127  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
1128 
1129  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
1130  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraNcnstCrsMat));
1131  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1132  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> ret =
1133  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat);
1134  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1135  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
1136  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1137  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1138  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1139  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1140  return xpMat;
1141 #else
1142  throw Xpetra::Exceptions::RuntimeError("Problem BBB. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1143 #endif
1144  }
1145 #endif
1146 
1147 #ifdef HAVE_XPETRA_EPETRA
1148  if (isEpetra(op)) {
1149  Teuchos::RCP<const Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*op);
1150  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1151  Teuchos::RCP<const Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_op, true);
1152  Teuchos::RCP<const Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(epetra_rowmat, true);
1153  Teuchos::RCP<Epetra_CrsMatrix> epetra_ncnstcrsmat = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(epetra_crsmat);
1154  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_ncnstcrsmat));
1155 
1156  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpetraCrsMat =
1157  Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>(epetra_ncnstcrsmat));
1158  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1159 
1160  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1161  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xEpetraCrsMat, true);
1162  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1163  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1164  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1165  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1166  return xpMat;
1167  }
1168 #endif
1169  return Teuchos::null;
1170  }
1171 
1172  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1173  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
1174 #ifdef HAVE_XPETRA_TPETRA
1175  if (isTpetra(op)) {
1176 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1177  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1178 
1179  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1180  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
1181  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1182  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
1183  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
1184 
1185  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
1186  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraCrsMat));
1187  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1188  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1189  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
1190  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1191  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1192  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1193  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1194  return xpMat;
1195 #else
1196  throw Xpetra::Exceptions::RuntimeError("Problem CCC. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1197 #endif
1198  }
1199 #endif
1200 
1201 #ifdef HAVE_XPETRA_EPETRA
1202  if (isEpetra(op)) {
1203  Teuchos::RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*op);
1204  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1205  Teuchos::RCP<Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetra_op, true);
1206  Teuchos::RCP<Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetra_rowmat, true);
1207 
1208  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpetraCrsMat =
1209  Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>(epetra_crsmat));
1210  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1211  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1212  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xEpetraCrsMat, true);
1213  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1214  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1215  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1216  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1217  return xpMat;
1218  }
1219 #endif
1220  return Teuchos::null;
1221  }
1222 
1223  static Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1224  toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1225  return toXpetraOperator(Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar>>(op));
1226 
1227  // #ifdef HAVE_XPETRA_TPETRA
1228  // if(isTpetra(op)) {
1229  // typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1230  // Teuchos::RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getConstTpetraOperator(op);
1231 
1232  // Teuchos::RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > nonConstTpetraOp =
1233  // Teuchos::rcp_const_cast<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
1234 
1235  // Teuchos::RCP<Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xTpetraOp =
1236  // Teuchos::rcp(new Xpetra::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node>(nonConstTpetraOp));
1237  // TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
1238 
1239  // Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xpOp =
1240  // Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xTpetraOp, true);
1241  // return xpOp;
1242  // }
1243  // #endif
1244 
1245  // #ifdef HAVE_XPETRA_EPETRA
1246  // if(isEpetra(op)) {
1247  // TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
1248  // }
1249  // #endif
1250  // return Teuchos::null;
1251  }
1252 
1253  static Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1254  toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
1255 #ifdef HAVE_XPETRA_TPETRA
1256  if (isTpetra(op)) {
1257  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1258  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
1259 
1260  Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
1261  Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
1262  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
1263 
1264  Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
1265  Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
1266  return xpOp;
1267  }
1268 #endif
1269 
1270 #ifdef HAVE_XPETRA_EPETRA
1271  if (isEpetra(op)) {
1272  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
1273  }
1274 #endif
1275  return Teuchos::null;
1276  }
1277 
1278  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1279  toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op) {
1280  using Teuchos::rcp_const_cast;
1281  using Teuchos::rcp_dynamic_cast;
1282 
1283  RCP<const Thyra::VectorBase<Scalar>> diag = op->getDiag();
1284 
1285  RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpDiag;
1286 #ifdef HAVE_XPETRA_TPETRA
1287  using thyTpV = Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
1288  using tV = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
1289  if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
1290  RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
1291  if (!tDiag.is_null())
1292  xpDiag = Xpetra::toXpetra(tDiag);
1293  }
1294 #endif
1295 #ifdef HAVE_XPETRA_EPETRA
1296  using ThyVSBase = Thyra::SpmdVectorSpaceBase<Scalar>;
1297  if (xpDiag.is_null()) {
1298  RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(op->range())->getComm());
1299  RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(op->range()), comm);
1300  if (!map.is_null()) {
1301  RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
1302  RCP<Epetra_Vector> nceDiag = rcp_const_cast<Epetra_Vector>(eDiag);
1303  RCP<Xpetra::EpetraVectorT<int, Node>> xpEpDiag = rcp(new Xpetra::EpetraVectorT<int, Node>(nceDiag));
1304  xpDiag = rcp_dynamic_cast<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpEpDiag, true);
1305  }
1306  }
1307 #endif
1308  TEUCHOS_ASSERT(!xpDiag.is_null());
1309  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
1310  return M;
1311  }
1312 
1313  static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1314  toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op) {
1315  return toXpetra(Teuchos::rcp_const_cast<Thyra::DiagonalLinearOpBase<Scalar>>(op));
1316  }
1317 
1318  static Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
1319  toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map) {
1320  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> thyraMap = Teuchos::null;
1321 
1322  // check whether map is of type BlockedMap
1323  RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
1324  if (bmap.is_null() == false) {
1325  Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>> vecSpaces(bmap->getNumMaps());
1326  for (size_t i = 0; i < bmap->getNumMaps(); i++) {
1327  // we need Thyra GIDs for all the submaps
1328  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> vs =
1329  Xpetra::ThyraUtils<Scalar, LO, GO, Node>::toThyra(bmap->getMap(i, true));
1330  vecSpaces[i] = vs;
1331  }
1332 
1333  thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
1334  return thyraMap;
1335  }
1336 
1337  // standard case
1338 #ifdef HAVE_XPETRA_TPETRA
1339  if (map->lib() == Xpetra::UseTpetra) {
1340 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1341  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1342  Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>> tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(map);
1343  if (tpetraMap == Teuchos::null)
1344  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
1345  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tpMap = tpetraMap->getTpetra_Map();
1346  RCP<Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
1347  thyraMap = thyraTpetraMap;
1348 #else
1349  throw Xpetra::Exceptions::RuntimeError("Problem DDD. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1350 #endif
1351  }
1352 #endif
1353 
1354 #ifdef HAVE_XPETRA_EPETRA
1355  if (map->lib() == Xpetra::UseEpetra) {
1356  Teuchos::RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>> epetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(map);
1357  if (epetraMap == Teuchos::null)
1358  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
1359  RCP<const Thyra::VectorSpaceBase<Scalar>> thyraEpetraMap = Thyra::create_VectorSpace(epetraMap->getEpetra_MapRCP());
1360  thyraMap = thyraEpetraMap;
1361  }
1362 #endif
1363 
1364  return thyraMap;
1365  }
1366 
1367  static Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
1368  toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
1369  // create Thyra MultiVector
1370 #ifdef HAVE_XPETRA_TPETRA
1371  if (vec->getMap()->lib() == Xpetra::UseTpetra) {
1372  auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
1373  RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpMV = Teuchos::rcp_dynamic_cast<const TpetraMultiVector>(vec)->getTpetra_MultiVector();
1374  auto thyDomMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(vec->getNumVectors(), vec->getMap()->getComm()));
1375  auto thyMV = rcp(new Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
1376  thyMV->initialize(thyTpMap, thyDomMap, tpMV);
1377  return thyMV;
1378  }
1379 #endif
1380 
1381 #ifdef HAVE_XPETRA_EPETRA
1382  if (vec->getMap()->lib() == Xpetra::UseEpetra) {
1383  auto thyEpMap = Thyra::create_VectorSpace(Teuchos::rcp_dynamic_cast<const EpetraMapT<GlobalOrdinal, Node>>(vec->getMap())->getEpetra_MapRCP());
1384  auto epMV = Teuchos::rcp_dynamic_cast<const EpetraMultiVectorT<GlobalOrdinal, Node>>(vec)->getEpetra_MultiVector();
1385  auto thyMV = Thyra::create_MultiVector(epMV, thyEpMap);
1386  return thyMV;
1387  }
1388 #endif
1389 
1390  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "MultiVector cannot be converted to Thyra.");
1391  }
1392 
1393  static Teuchos::RCP<const Thyra::VectorBase<Scalar>>
1394  toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
1395  // create Thyra Vector
1396 #ifdef HAVE_XPETRA_TPETRA
1397  if (vec->getMap()->lib() == Xpetra::UseTpetra) {
1398  auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
1399  RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpVec = Teuchos::rcp_dynamic_cast<const TpetraVector>(vec)->getTpetra_Vector();
1400  auto thyVec = rcp(new Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
1401  thyVec->initialize(thyTpMap, tpVec);
1402  return thyVec;
1403  }
1404 #endif
1405 
1406 #ifdef HAVE_XPETRA_EPETRA
1407  if (vec->getMap()->lib() == Xpetra::UseEpetra) {
1408  auto thyEpMap = Thyra::create_VectorSpace(Teuchos::rcp_dynamic_cast<const EpetraMapT<GlobalOrdinal, Node>>(vec->getMap())->getEpetra_MapRCP());
1409  auto epVec = rcp(Teuchos::rcp_dynamic_cast<const EpetraVectorT<GlobalOrdinal, Node>>(vec)->getEpetra_Vector(), false);
1410  auto thyVec = Thyra::create_Vector(epVec, thyEpMap);
1411  return thyVec;
1412  }
1413 #endif
1414 
1415  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Vector cannot be converted to Thyra.");
1416  }
1417 
1418  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) {
1419  using Teuchos::as;
1420  using Teuchos::RCP;
1421  using Teuchos::rcp_dynamic_cast;
1422  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1423  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
1424  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1425  // typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
1426  // typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
1427  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1428 
1429  // copy data from tY_inout to Y_inout
1430  RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
1431  if (prodTarget != Teuchos::null) {
1432  RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
1433  if (bSourceVec.is_null() == true) {
1434  // SPECIAL CASE: target vector is product vector:
1435  // update Thyra product multi vector with data from a merged Xpetra multi vector
1436 
1437  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
1438  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.");
1439 
1440  for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1441  // access Xpetra data
1442  RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
1443 
1444  // access Thyra data
1445  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
1446  RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
1447  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
1448  const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
1449  const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim());
1450  RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
1451  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
1452 
1453  // loop over all vectors in multivector
1454  for (size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
1455  Teuchos::ArrayRCP<const Scalar> xpData = xpSubBlock->getData(j); // access const data from Xpetra object
1456 
1457  // loop over all local rows
1458  for (LocalOrdinal i = 0; i < localSubDim; ++i) {
1459  (*thyData)(i, j) = xpData[i];
1460  }
1461  }
1462  }
1463  } else {
1464  // source vector is a blocked multivector
1465  // TODO test me
1466  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.");
1467 
1468  for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1469  // access Thyra data
1470  RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
1471 
1472  Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
1473 
1474  // access Thyra data
1475  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
1476  Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
1477  }
1478  }
1479  } else {
1480  // STANDARD case:
1481  // update Thyra::MultiVector with data from an Xpetra::MultiVector
1482 
1483  // access Thyra data
1484  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
1485  TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
1486  const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
1487  const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim());
1488  RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
1489  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
1490 
1491  // loop over all vectors in multivector
1492  for (size_t j = 0; j < source->getNumVectors(); ++j) {
1493  Teuchos::ArrayRCP<const Scalar> xpData = source->getData(j); // access const data from Xpetra object
1494  // loop over all local rows
1495  for (LocalOrdinal i = 0; i < localSubDim; ++i) {
1496  (*thyData)(i, j) = xpData[i];
1497  }
1498  }
1499  }
1500  }
1501 
1502  static Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
1503  toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
1504  // create a Thyra operator from Xpetra::CrsMatrix
1505  Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
1506 
1507 #ifdef HAVE_XPETRA_TPETRA
1508  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
1509  if (tpetraMat != Teuchos::null) {
1510 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1511  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1512 
1513  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
1514  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
1515  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
1516 
1517  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
1518  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
1519 
1520  thyraOp = Thyra::createConstLinearOp(tpOperator);
1521  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
1522 #else
1523  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1524 #endif
1525  }
1526 #endif
1527 
1528 #ifdef HAVE_XPETRA_EPETRA
1529  Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
1530  if (epetraMat != Teuchos::null) {
1531  Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
1532  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpCrsMat));
1533  Teuchos::RCP<const Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrix();
1534  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
1535 
1536  Teuchos::RCP<const Thyra::EpetraLinearOp> thyraEpOp = Thyra::epetraLinearOp(epCrsMat, "op");
1537  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
1538  thyraOp = thyraEpOp;
1539  }
1540 #endif
1541  return thyraOp;
1542  }
1543 
1544  static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
1545  toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
1546  // create a Thyra operator from Xpetra::CrsMatrix
1547  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
1548 
1549 #ifdef HAVE_XPETRA_TPETRA
1550  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
1551  if (tpetraMat != Teuchos::null) {
1552 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1553  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1554 
1555  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
1556  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
1557  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
1558 
1559  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
1560  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
1561 
1562  thyraOp = Thyra::createLinearOp(tpOperator);
1563  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
1564 #else
1565  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1566 #endif
1567  }
1568 #endif
1569 
1570 #ifdef HAVE_XPETRA_EPETRA
1571  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
1572  if (epetraMat != Teuchos::null) {
1573  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat, true);
1574  Teuchos::RCP<Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrixNonConst();
1575  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
1576 
1577  Teuchos::RCP<Thyra::EpetraLinearOp> thyraEpOp = Thyra::nonconstEpetraLinearOp(epCrsMat, "op");
1578  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
1579  thyraOp = thyraEpOp;
1580  }
1581 #endif
1582  return thyraOp;
1583  }
1584 
1585  static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
1586  toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, int, EpetraNode>>& mat);
1587 
1588 }; // specialization on SC=double, LO=GO=int and NO=EpetraNode
1589 #endif // #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
1590 
1591 #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
1592 template <>
1593 class ThyraUtils<double, int, long long, EpetraNode> {
1594  public:
1595  typedef double Scalar;
1596  typedef int LocalOrdinal;
1597  typedef long long GlobalOrdinal;
1598  typedef EpetraNode Node;
1599 
1600  private:
1601 #undef XPETRA_THYRAUTILS_SHORT
1602 #include "Xpetra_UseShortNames.hpp"
1603 
1604  public:
1605  static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>
1606  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) {
1607  Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map = ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(vectorSpace, comm);
1608 
1609  if (stridedBlockId == -1) {
1610  TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo.size() != 0);
1611  } else {
1612  TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo[stridedBlockId] != 0);
1613  }
1614 
1615  Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>> ret = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(map, stridingInfo, stridedBlockId, offset);
1616  return ret;
1617  }
1618 
1619  static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
1620  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
1621  using Teuchos::as;
1622  using Teuchos::RCP;
1623  using Teuchos::rcp_dynamic_cast;
1624  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1625  typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
1626  typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyUtils;
1627 
1628  RCP<const ThyProdVecSpaceBase> prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase>(vectorSpace);
1629  if (prodVectorSpace != Teuchos::null) {
1630  // SPECIAL CASE: product Vector space
1631  // collect all submaps to store them in a hierarchical BlockedMap object
1632  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks() == 0, std::logic_error, "Found a product vector space with zero blocks.");
1633  std::vector<RCP<const Map>> mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
1634  std::vector<RCP<const Map>> mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
1635  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
1636  RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
1637  // can be of type Map or BlockedMap (containing Thyra GIDs)
1638  mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
1639  }
1640 
1641  // get offsets for submap GIDs
1642  // we need that for the full map (Xpetra GIDs)
1643  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(), 0);
1644  for (int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
1645  gidOffsets[i] = mapsThyra[i - 1]->getMaxAllGlobalIndex() + gidOffsets[i - 1] + 1;
1646  }
1647 
1648  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
1649  RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
1650  // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
1651  mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
1652  }
1653 
1654  Teuchos::RCP<Map> resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapsXpetra, mapsThyra));
1655  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
1656  return resultMap;
1657  } else {
1658  // STANDARD CASE: no product map
1659  // Epetra/Tpetra specific code to access the underlying map data
1660 
1661  // check whether we have a Tpetra based Thyra operator
1662  bool bIsTpetra = false;
1663 #ifdef HAVE_XPETRA_TPETRA
1664 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1665  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1666  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vectorSpace);
1667  bIsTpetra = Teuchos::is_null(tpetra_vsc) ? false : true;
1668 #endif
1669 #endif
1670 
1671  // check whether we have an Epetra based Thyra operator
1672  bool bIsEpetra = !bIsTpetra; // note: this is a little bit fragile!
1673 
1674 #ifdef HAVE_XPETRA_TPETRA
1675  if (bIsTpetra) {
1676 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1677  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1678  typedef Thyra::VectorBase<Scalar> ThyVecBase;
1679  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> TpMap;
1680  typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpVector;
1681  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1682  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
1683  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
1684  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
1685  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
1686  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
1687  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
1688 
1689  RCP<Map> rgXpetraMap = Xpetra::toXpetraNonConst(rgTpetraMap);
1690  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgXpetraMap));
1691  return rgXpetraMap;
1692 #else
1693  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1694 #endif
1695  }
1696 #endif
1697 
1698 #ifdef HAVE_XPETRA_EPETRA
1699  if (bIsEpetra) {
1700  // RCP<const Epetra_Map> epMap = Teuchos::null;
1701  RCP<const Epetra_Map>
1702  epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map>>(vectorSpace, "epetra_map");
1703  if (!Teuchos::is_null(epetra_map)) {
1704  Teuchos::RCP<Map> rgXpetraMap = Teuchos::rcp(new Xpetra::EpetraMapT<GlobalOrdinal, Node>(epetra_map));
1705  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgXpetraMap));
1706  return rgXpetraMap;
1707  } else {
1708  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "No Epetra_Map data found in Thyra::VectorSpace.");
1709  }
1710  }
1711 #endif
1712  } // end standard case (no product map)
1713  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Cannot transform Thyra::VectorSpace to Xpetra::Map.");
1714  // return Teuchos::null; // unreachable
1715  }
1716 
1717  // const version
1718  static Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1719  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
1720  using Teuchos::as;
1721  using Teuchos::RCP;
1722  using Teuchos::rcp_dynamic_cast;
1723  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1724  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1725  typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyUtils;
1726 
1727  // return value
1728  RCP<MultiVector> xpMultVec = Teuchos::null;
1729 
1730  // check whether v is a product multi vector
1731  Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase>(v);
1732  if (thyProdVec != Teuchos::null) {
1733  // SPECIAL CASE: create a nested BlockedMultiVector
1734  // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
1735  RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
1736  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
1737 
1738  // create new Xpetra::BlockedMultiVector
1739  xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
1740 
1741  RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec);
1742  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpBlockedMultVec));
1743 
1744  // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
1745  for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b) {
1746  RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
1747  // xpBlockMV can be of type MultiVector or BlockedMultiVector
1748  RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); // recursive call
1749  xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
1750  }
1751 
1752  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1753  return xpMultVec;
1754  } else {
1755  // STANDARD CASE: no product vector
1756  // Epetra/Tpetra specific code to access the underlying map data
1757  bool bIsTpetra = false;
1758 #ifdef HAVE_XPETRA_TPETRA
1759 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1760  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1761 
1762  // typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
1763  // typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
1764  typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
1765  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> ConverterT;
1766  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpMultVec;
1767  typedef Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpTpMultVec;
1768  typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
1769 
1770  RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
1771  RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
1772  if (thyraTpetraMultiVector != Teuchos::null) {
1773  bIsTpetra = true;
1774  const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
1775  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
1776  RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
1777  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
1778  xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
1779  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1780  return xpMultVec;
1781  }
1782 #endif
1783 #endif
1784 
1785 #ifdef HAVE_XPETRA_EPETRA
1786  if (bIsTpetra == false) {
1787  // no product vector
1788  Teuchos::RCP<Map> map = ThyUtils::toXpetra(v->range(), comm);
1789  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(map));
1790  RCP<Xpetra::EpetraMapT<GlobalOrdinal, Node>> xeMap = rcp_dynamic_cast<Xpetra::EpetraMapT<GlobalOrdinal, Node>>(map, true);
1791  RCP<const Epetra_Map> eMap = xeMap->getEpetra_MapRCP();
1792  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(eMap));
1793  Teuchos::RCP<const Epetra_MultiVector> epMultVec = Thyra::get_Epetra_MultiVector(*eMap, v);
1794  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epMultVec));
1795  RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<Epetra_MultiVector>(epMultVec);
1796  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epNonConstMultVec));
1797  xpMultVec = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GlobalOrdinal, Node>(epNonConstMultVec));
1798  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1799  return xpMultVec;
1800  }
1801 #endif
1802  } // end standard case
1803  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Cannot transform Thyra::MultiVector to Xpetra::MultiVector.");
1804  // return Teuchos::null; // unreachable
1805  }
1806 
1807  // non-const version
1808  static Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1809  toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
1810  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> cv =
1811  Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar>>(v);
1812  Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> r =
1813  toXpetra(cv, comm);
1814  return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(r);
1815  }
1816 
1817  static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1818  // check whether we have a Tpetra based Thyra operator
1819  bool bIsTpetra = false;
1820 #ifdef HAVE_XPETRA_TPETRA
1821 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1822  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1823 
1824  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(op);
1825  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
1826 
1827  // for debugging purposes: find out why dynamic cast failed
1828  if (!bIsTpetra &&
1829 #ifdef HAVE_XPETRA_EPETRA
1830  Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
1831 #endif
1832  Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar>>(op) == Teuchos::null) {
1833  // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
1834  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
1835  if (Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
1836  std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
1837  std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
1838  std::cout << " properly set!" << std::endl;
1839  std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
1840  }
1841  }
1842 #endif
1843 #endif
1844 
1845 #if 0
1846  // Check whether it is a blocked operator.
1847  // If yes, grab the (0,0) block and check the underlying linear algebra
1848  // Note: we expect that the (0,0) block exists!
1849  if(bIsTpetra == false) {
1850  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1851  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
1852  if(ThyBlockedOp != Teuchos::null) {
1853  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1854  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1855  ThyBlockedOp->getBlock(0,0);
1856  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1857  bIsTpetra = isTpetra(b00);
1858  }
1859  }
1860 #endif
1861 
1862  return bIsTpetra;
1863  }
1864 
1865  static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1866  // check whether we have an Epetra based Thyra operator
1867  bool bIsEpetra = false;
1868 
1869 #ifdef HAVE_XPETRA_EPETRA
1870  Teuchos::RCP<const Thyra::EpetraLinearOp> epetra_op = Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op, false);
1871  bIsEpetra = Teuchos::is_null(epetra_op) ? false : true;
1872 #endif
1873 
1874 #if 0
1875  // Check whether it is a blocked operator.
1876  // If yes, grab the (0,0) block and check the underlying linear algebra
1877  // Note: we expect that the (0,0) block exists!
1878  if(bIsEpetra == false) {
1879  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1880  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op,false);
1881  if(ThyBlockedOp != Teuchos::null) {
1882  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1883  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1884  ThyBlockedOp->getBlock(0,0);
1885  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1886  bIsEpetra = isEpetra(b00);
1887  }
1888  }
1889 #endif
1890 
1891  return bIsEpetra;
1892  }
1893 
1894  static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1895  // Check whether it is a blocked operator.
1896  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar>> ThyBlockedOp =
1897  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar>>(op);
1898  if (ThyBlockedOp != Teuchos::null) {
1899  return true;
1900  }
1901  return false;
1902  }
1903 
1904  static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1905  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
1906 #ifdef HAVE_XPETRA_TPETRA
1907  if (isTpetra(op)) {
1908 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1909  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1910 
1911  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1912  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
1913  // we should also add support for the const versions!
1914  // getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
1915  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1916  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
1917  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
1918  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat);
1919  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
1920 
1921  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
1922  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraNcnstCrsMat));
1923  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1924 
1925  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1926  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
1927  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1928  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1929  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1930  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1931  return xpMat;
1932 #else
1933  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1934 #endif
1935  }
1936 #endif
1937 
1938 #ifdef HAVE_XPETRA_EPETRA
1939  if (isEpetra(op)) {
1940  Teuchos::RCP<const Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*op);
1941  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1942  Teuchos::RCP<const Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_op, true);
1943  Teuchos::RCP<const Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(epetra_rowmat, true);
1944  Teuchos::RCP<Epetra_CrsMatrix> epetra_ncnstcrsmat = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(epetra_crsmat);
1945  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_ncnstcrsmat));
1946 
1947  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpetraCrsMat =
1948  Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>(epetra_ncnstcrsmat));
1949  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1950 
1951  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1952  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xEpetraCrsMat, true);
1953  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1954  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1955  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1956  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1957  return xpMat;
1958  }
1959 #endif
1960  return Teuchos::null;
1961  }
1962 
1963  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
1964  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
1965 #ifdef HAVE_XPETRA_TPETRA
1966  if (isTpetra(op)) {
1967 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1968  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1969 
1970  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
1971  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
1972  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1973  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
1974  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
1975 
1976  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
1977  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraCrsMat));
1978  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1979 
1980  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
1981  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
1982  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
1983  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
1984  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
1985  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
1986  return xpMat;
1987 #else
1988  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1989 #endif
1990  }
1991 #endif
1992 
1993 #ifdef HAVE_XPETRA_EPETRA
1994  if (isEpetra(op)) {
1995  Teuchos::RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*op);
1996  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1997  Teuchos::RCP<Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetra_op, true);
1998  Teuchos::RCP<Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetra_rowmat, true);
1999 
2000  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpetraCrsMat =
2001  Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>(epetra_crsmat));
2002  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
2003 
2004  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
2005  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xEpetraCrsMat, true);
2006  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
2007  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
2008  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
2009  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
2010  return xpMat;
2011  }
2012 #endif
2013  return Teuchos::null;
2014  }
2015 
2016  static Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
2017  toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
2018 #ifdef HAVE_XPETRA_TPETRA
2019  if (isTpetra(op)) {
2020  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
2021  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
2022 
2023  Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
2024  Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
2025  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
2026 
2027  Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
2028  Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
2029  return xpOp;
2030  }
2031 #endif
2032 
2033 #ifdef HAVE_XPETRA_EPETRA
2034  if (isEpetra(op)) {
2035  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
2036  }
2037 #endif
2038  return Teuchos::null;
2039  }
2040 
2041  static Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
2042  toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
2043 #ifdef HAVE_XPETRA_TPETRA
2044  if (isTpetra(op)) {
2045  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
2046  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
2047 
2048  Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
2049  Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
2050  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
2051 
2052  Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
2053  Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
2054  return xpOp;
2055  }
2056 #endif
2057 
2058 #ifdef HAVE_XPETRA_EPETRA
2059  if (isEpetra(op)) {
2060  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
2061  }
2062 #endif
2063  return Teuchos::null;
2064  }
2065 
2066  static Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
2067  toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op) {
2068  using Teuchos::rcp_const_cast;
2069  using Teuchos::rcp_dynamic_cast;
2070  using ThyVSBase = Thyra::SpmdVectorSpaceBase<Scalar>;
2071  using thyTpV = Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
2072  using tV = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
2073 
2074  RCP<const Thyra::VectorBase<Scalar>> diag = op->getDiag();
2075 
2076  RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpDiag;
2077 #ifdef HAVE_XPETRA_TPETRA
2078  if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
2079  RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
2080  if (!tDiag.is_null())
2081  xpDiag = Xpetra::toXpetra(tDiag);
2082  }
2083 #endif
2084 #ifdef HAVE_XPETRA_EPETRA
2085  if (xpDiag.is_null()) {
2086  RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(op->range())->getComm());
2087  RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(op->range()), comm);
2088  if (!map.is_null()) {
2089  RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
2090  RCP<Epetra_Vector> nceDiag = rcp_const_cast<Epetra_Vector>(eDiag);
2091  RCP<Xpetra::EpetraVectorT<int, Node>> xpEpDiag = rcp(new Xpetra::EpetraVectorT<int, Node>(nceDiag));
2092  xpDiag = rcp_dynamic_cast<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpEpDiag, true);
2093  }
2094  }
2095 #endif
2096  TEUCHOS_ASSERT(!xpDiag.is_null());
2097  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
2098  return M;
2099  }
2100 
2101  static Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
2102  toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op) {
2103  return toXpetra(Teuchos::rcp_const_cast<Thyra::DiagonalLinearOpBase<Scalar>>(op));
2104  }
2105 
2106  static Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
2107  toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map) {
2108  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> thyraMap = Teuchos::null;
2109 
2110  // check whether map is of type BlockedMap
2111  RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
2112  if (bmap.is_null() == false) {
2113  Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>> vecSpaces(bmap->getNumMaps());
2114  for (size_t i = 0; i < bmap->getNumMaps(); i++) {
2115  // we need Thyra GIDs for all the submaps
2116  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> vs =
2117  Xpetra::ThyraUtils<Scalar, LO, GO, Node>::toThyra(bmap->getMap(i, true));
2118  vecSpaces[i] = vs;
2119  }
2120 
2121  thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
2122  return thyraMap;
2123  }
2124 
2125  // standard case
2126 #ifdef HAVE_XPETRA_TPETRA
2127  if (map->lib() == Xpetra::UseTpetra) {
2128 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
2129  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
2130  Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>> tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(map);
2131  if (tpetraMap == Teuchos::null)
2132  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
2133  RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> tpMap = tpetraMap->getTpetra_Map();
2134  RCP<Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
2135  thyraMap = thyraTpetraMap;
2136 #else
2137  throw Xpetra::Exceptions::RuntimeError("Problem DDD. Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
2138 #endif
2139  }
2140 #endif
2141 
2142 #ifdef HAVE_XPETRA_EPETRA
2143  if (map->lib() == Xpetra::UseEpetra) {
2144  Teuchos::RCP<const Xpetra::EpetraMapT<GlobalOrdinal, Node>> epetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal, Node>>(map);
2145  if (epetraMap == Teuchos::null)
2146  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
2147  RCP<const Thyra::VectorSpaceBase<Scalar>> thyraEpetraMap = Thyra::create_VectorSpace(epetraMap->getEpetra_MapRCP());
2148  thyraMap = thyraEpetraMap;
2149  }
2150 #endif
2151 
2152  return thyraMap;
2153  }
2154 
2155  static Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
2156  toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
2157  // create Thyra MultiVector
2158 #ifdef HAVE_XPETRA_TPETRA
2159  if (vec->getMap()->lib() == Xpetra::UseTpetra) {
2160  auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
2161  RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpMV = Teuchos::rcp_dynamic_cast<const TpetraMultiVector>(vec)->getTpetra_MultiVector();
2162  auto thyDomMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(vec->getNumVectors(), vec->getMap()->getComm()));
2163  auto thyMV = rcp(new Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
2164  thyMV->initialize(thyTpMap, thyDomMap, tpMV);
2165  return thyMV;
2166  }
2167 #endif
2168 
2169 #ifdef HAVE_XPETRA_EPETRA
2170  if (vec->getMap()->lib() == Xpetra::UseEpetra) {
2171  auto thyEpMap = Thyra::create_VectorSpace(Teuchos::rcp_dynamic_cast<const EpetraMapT<GlobalOrdinal, Node>>(vec->getMap())->getEpetra_MapRCP());
2172  auto epMV = Teuchos::rcp_dynamic_cast<const EpetraMultiVectorT<GlobalOrdinal, Node>>(vec)->getEpetra_MultiVector();
2173  auto thyMV = Thyra::create_MultiVector(epMV, thyEpMap);
2174  return thyMV;
2175  }
2176 #endif
2177 
2178  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "MultiVector cannot be converted to Thyra.");
2179  }
2180 
2181  static Teuchos::RCP<const Thyra::VectorBase<Scalar>>
2182  toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
2183  // create Thyra Vector
2184 #ifdef HAVE_XPETRA_TPETRA
2185  if (vec->getMap()->lib() == Xpetra::UseTpetra) {
2186  auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
2187  RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpVec = Teuchos::rcp_dynamic_cast<const TpetraVector>(vec)->getTpetra_Vector();
2188  auto thyVec = rcp(new Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
2189  thyVec->initialize(thyTpMap, tpVec);
2190  return thyVec;
2191  }
2192 #endif
2193 
2194 #ifdef HAVE_XPETRA_EPETRA
2195  if (vec->getMap()->lib() == Xpetra::UseEpetra) {
2196  auto thyEpMap = Thyra::create_VectorSpace(Teuchos::rcp_dynamic_cast<const EpetraMapT<GlobalOrdinal, Node>>(vec->getMap())->getEpetra_MapRCP());
2197  auto epVec = rcp(Teuchos::rcp_dynamic_cast<const EpetraVectorT<GlobalOrdinal, Node>>(vec)->getEpetra_Vector(), false);
2198  auto thyVec = Thyra::create_Vector(epVec, thyEpMap);
2199  return thyVec;
2200  }
2201 #endif
2202 
2203  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Vector cannot be converted to Thyra.");
2204  }
2205 
2206  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) {
2207  using Teuchos::as;
2208  using Teuchos::RCP;
2209  using Teuchos::rcp_dynamic_cast;
2210  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
2211  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
2212  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
2213  // typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
2214  // typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
2215  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
2216 
2217  // copy data from tY_inout to Y_inout
2218  RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
2219  if (prodTarget != Teuchos::null) {
2220  RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
2221  if (bSourceVec.is_null() == true) {
2222  // SPECIAL CASE: target vector is product vector:
2223  // update Thyra product multi vector with data from a merged Xpetra multi vector
2224 
2225  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
2226  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.");
2227 
2228  for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
2229  // access Xpetra data
2230  RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
2231 
2232  // access Thyra data
2233  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
2234  RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
2235  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
2236  const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
2237  const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim());
2238  RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
2239  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
2240 
2241  // loop over all vectors in multivector
2242  for (size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
2243  Teuchos::ArrayRCP<const Scalar> xpData = xpSubBlock->getData(j); // access const data from Xpetra object
2244 
2245  // loop over all local rows
2246  for (LocalOrdinal i = 0; i < localSubDim; ++i) {
2247  (*thyData)(i, j) = xpData[i];
2248  }
2249  }
2250  }
2251  } else {
2252  // source vector is a blocked multivector
2253  // TODO test me
2254  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.");
2255 
2256  for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
2257  // access Thyra data
2258  RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
2259 
2260  Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
2261 
2262  // access Thyra data
2263  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
2264  Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
2265  }
2266  }
2267  } else {
2268  // STANDARD case:
2269  // update Thyra::MultiVector with data from an Xpetra::MultiVector
2270 
2271  // access Thyra data
2272  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
2273  TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
2274  const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
2275  const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim());
2276  RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
2277  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
2278 
2279  // loop over all vectors in multivector
2280  for (size_t j = 0; j < source->getNumVectors(); ++j) {
2281  Teuchos::ArrayRCP<const Scalar> xpData = source->getData(j); // access const data from Xpetra object
2282  // loop over all local rows
2283  for (LocalOrdinal i = 0; i < localSubDim; ++i) {
2284  (*thyData)(i, j) = xpData[i];
2285  }
2286  }
2287  }
2288  }
2289 
2290  static Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
2291  toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
2292  // create a Thyra operator from Xpetra::CrsMatrix
2293  Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
2294 
2295 #ifdef HAVE_XPETRA_TPETRA
2296  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
2297  if (tpetraMat != Teuchos::null) {
2298 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
2299  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
2300 
2301  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
2302  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
2303  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
2304 
2305  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
2306  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
2307 
2308  thyraOp = Thyra::createConstLinearOp(tpOperator);
2309  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
2310 #else
2311  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
2312 #endif
2313  }
2314 #endif
2315 
2316 #ifdef HAVE_XPETRA_EPETRA
2317  Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
2318  if (epetraMat != Teuchos::null) {
2319  Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat, true);
2320  Teuchos::RCP<const Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrix();
2321  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
2322 
2323  Teuchos::RCP<const Thyra::EpetraLinearOp> thyraEpOp = Thyra::epetraLinearOp(epCrsMat, "op");
2324  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
2325  thyraOp = thyraEpOp;
2326  }
2327 #endif
2328  return thyraOp;
2329  }
2330 
2331  static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
2332  toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
2333  // create a Thyra operator from Xpetra::CrsMatrix
2334  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
2335 
2336 #ifdef HAVE_XPETRA_TPETRA
2337  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
2338  if (tpetraMat != Teuchos::null) {
2339 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
2340  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
2341 
2342  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
2343  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
2344  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
2345 
2346  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
2347  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
2348 
2349  thyraOp = Thyra::createLinearOp(tpOperator);
2350  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
2351 #else
2352  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
2353 #endif
2354  }
2355 #endif
2356 
2357 #ifdef HAVE_XPETRA_EPETRA
2358  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat);
2359  if (epetraMat != Teuchos::null) {
2360  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> xEpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(mat, true);
2361  Teuchos::RCP<Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrixNonConst();
2362  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
2363 
2364  Teuchos::RCP<Thyra::EpetraLinearOp> thyraEpOp = Thyra::nonconstEpetraLinearOp(epCrsMat, "op");
2365  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
2366  thyraOp = thyraEpOp;
2367  }
2368 #endif
2369  return thyraOp;
2370  }
2371 
2372  static Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
2373  toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, long long, EpetraNode>>& mat);
2374 
2375 }; // specialization on SC=double, LO=GO=int and NO=EpetraNode
2376 #endif // XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
2377 
2378 #endif // HAVE_XPETRA_EPETRA
2379 
2380 } // end namespace Xpetra
2381 
2382 #define XPETRA_THYRAUTILS_SHORT
2383 #endif // HAVE_XPETRA_THYRA
2384 
2385 #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:60
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrixNonConst
RCP< Epetra_CrsMatrix > getEpetra_CrsMatrixNonConst() const
Definition: Xpetra_EpetraCrsMatrix.hpp:242
Xpetra_UseShortNames.hpp
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory.hpp:235
Xpetra_MapUtils.hpp
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:101
Xpetra_Matrix.hpp
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:228
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrix
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
Definition: Xpetra_EpetraCrsMatrix.hpp:241
Xpetra_EpetraCrsMatrix.hpp
Xpetra_EpetraMap.hpp
EpetraNode
Tpetra::KokkosCompat::KokkosSerialWrapperNode EpetraNode
Definition: Xpetra_ETI_3arg.hpp:22
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
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:647
Xpetra::EpetraMapT::getEpetra_MapRCP
const RCP< const Epetra_Map > & getEpetra_MapRCP() const
Get the underlying Epetra map.
Definition: Xpetra_EpetraMap.hpp:246
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:651
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_fwd.hpp:51
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:61
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_fwd.hpp:51
Xpetra::toXpetraNonConst
const RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > toXpetraNonConst(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &map)
Definition: Xpetra_TpetraMap_decl.hpp:264
Xpetra::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Definition: Xpetra_EpetraCrsGraph.cpp:65