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