Xpetra  Version of the Day
Xpetra_ThyraUtils.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef XPETRA_THYRAUTILS_HPP
48 #define XPETRA_THYRAUTILS_HPP
49 
50 #include "Xpetra_ConfigDefs.hpp"
51 #ifdef HAVE_XPETRA_THYRA
52 
53 #include <typeinfo>
54 
55 #ifdef HAVE_XPETRA_TPETRA
56 #include "Tpetra_ConfigDefs.hpp"
57 #endif
58 
59 #ifdef HAVE_XPETRA_EPETRA
60 #include "Epetra_config.h"
61 #include "Epetra_CombineMode.h"
62 #endif
63 
64 #include "Xpetra_Map.hpp"
65 #include "Xpetra_BlockedMap.hpp"
66 #include "Xpetra_BlockedMultiVector.hpp"
67 #include "Xpetra_MapUtils.hpp"
68 #include "Xpetra_StridedMap.hpp"
69 #include "Xpetra_StridedMapFactory.hpp"
70 #include "Xpetra_MapExtractor.hpp"
71 #include "Xpetra_Matrix.hpp"
72 #include "Xpetra_CrsMatrixWrap.hpp"
73 #include "Xpetra_MultiVectorFactory.hpp"
74 
75 #include <Thyra_VectorSpaceBase.hpp>
76 #include <Thyra_SpmdVectorSpaceBase.hpp>
77 #include <Thyra_ProductVectorSpaceBase.hpp>
78 #include <Thyra_ProductMultiVectorBase.hpp>
79 #include <Thyra_VectorSpaceBase.hpp>
80 #include <Thyra_DefaultProductVectorSpace.hpp>
81 #include <Thyra_DefaultBlockedLinearOp.hpp>
82 #include <Thyra_LinearOpBase.hpp>
83 #include <Thyra_DetachedMultiVectorView.hpp>
84 #include <Thyra_MultiVectorStdOps.hpp>
85 
86 #ifdef HAVE_XPETRA_TPETRA
87 #include <Thyra_TpetraThyraWrappers.hpp>
88 #include <Thyra_TpetraVector.hpp>
89 #include <Thyra_TpetraVectorSpace.hpp>
90 #include <Tpetra_Map.hpp>
91 #include <Tpetra_Vector.hpp>
92 #include <Tpetra_CrsMatrix.hpp>
93 #include <Xpetra_TpetraMap.hpp>
94 #include <Xpetra_TpetraCrsMatrix.hpp>
95 #endif
96 #ifdef HAVE_XPETRA_EPETRA
97 #include <Thyra_EpetraLinearOp.hpp>
98 #include <Thyra_EpetraThyraWrappers.hpp>
99 #include <Thyra_SpmdVectorBase.hpp>
100 #include <Thyra_get_Epetra_Operator.hpp>
101 #include <Epetra_Map.h>
102 #include <Epetra_Vector.h>
103 #include <Epetra_CrsMatrix.h>
104 #include <Xpetra_EpetraMap.hpp>
105 #include <Xpetra_EpetraCrsMatrix.hpp>
106 #endif
107 
108 namespace Xpetra {
109 
110 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> class BlockedCrsMatrix;
111 
112 template <class Scalar,
113 class LocalOrdinal = int,
114 class GlobalOrdinal = LocalOrdinal,
115 class Node = KokkosClassic::DefaultNode::DefaultNodeType>
116 class ThyraUtils {
117 
118 private:
119 #undef XPETRA_THYRAUTILS_SHORT
120 #include "Xpetra_UseShortNames.hpp"
121 
122 public:
123  static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >
124  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) {
125 
126  Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map = toXpetra(vectorSpace);
127 
128  if(stridedBlockId == -1) {
129  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo.size() != 0);
130  }
131  else {
132  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo[stridedBlockId] != 0);
133  }
134 
135  Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> > ret = Xpetra::StridedMapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(map, stridingInfo, stridedBlockId, offset);
136  return ret;
137  }
138 
139  static Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
140  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
141  using Teuchos::RCP;
142  using Teuchos::rcp_dynamic_cast;
143  using Teuchos::as;
144  using ThyVecSpaceBase = Thyra::VectorSpaceBase<Scalar>;
145  using ThyProdVecSpaceBase = Thyra::ProductVectorSpaceBase<Scalar>;
146  using ThyUtils = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
147 
148 
149  RCP<Map> resultMap = Teuchos::null;
150  RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase >(vectorSpace);
151  if(prodVectorSpace != Teuchos::null) {
152  // SPECIAL CASE: product Vector space
153  // collect all submaps to store them in a hierarchical BlockedMap object
154  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error, "Found a product vector space with zero blocks.");
155  std::vector<RCP<const Map> > mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
156  std::vector<RCP<const Map> > mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
157  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
158  RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
159  // can be of type Map or BlockedMap (containing Thyra GIDs)
160  mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
161  }
162 
163  // get offsets for submap GIDs
164  // we need that for the full map (Xpetra GIDs)
165  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
166  for(int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
167  gidOffsets[i] = mapsThyra[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
168  }
169 
170  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
171  RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
172  // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
173  mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
174  }
175 
176  resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(mapsXpetra, mapsThyra));
177  } else {
178 #ifdef HAVE_XPETRA_TPETRA
179  // STANDARD CASE: no product map
180  // check whether we have a Tpetra based Thyra operator
181  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vectorSpace);
182  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.");
183 
184  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
185  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
186  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
187  typedef Thyra::VectorBase<Scalar> ThyVecBase;
188  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
189  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
190  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
191  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
192  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
193  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
194  resultMap = Xpetra::toXpetraNonConst(rgTpetraMap);
195 #else
196  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.");
197 #endif
198  } // end standard case (no product map)
199  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
200  return resultMap;
201  }
202 
203  // const version
204  static Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
205  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
206  using Teuchos::RCP;
207  using Teuchos::rcp_dynamic_cast;
208  using Teuchos::as;
209 
210  using ThyMultVecBase = Thyra::MultiVectorBase<Scalar>;
211  using ThyProdMultVecBase = Thyra::ProductMultiVectorBase<Scalar>;
212  using ThyUtils = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
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 
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::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
787 
788  RCP<Map> resultMap = Teuchos::null;
789 
790  RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase >(vectorSpace);
791  if(prodVectorSpace != Teuchos::null) {
792  // SPECIAL CASE: product Vector space
793  // collect all submaps to store them in a hierarchical BlockedMap object
794  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error, "Found a product vector space with zero blocks.");
795  std::vector<RCP<const Map> > mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
796  std::vector<RCP<const Map> > mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
797  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
798  RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
799  // can be of type Map or BlockedMap (containing Thyra GIDs)
800  mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
801  }
802 
803  // get offsets for submap GIDs
804  // we need that for the full map (Xpetra GIDs)
805  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
806  for(int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
807  gidOffsets[i] = mapsThyra[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
808  }
809 
810  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
811  RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
812  // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
813  mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
814  }
815 
816  resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(mapsXpetra, mapsThyra));
817  } else {
818  // STANDARD CASE: no product map
819  // Epetra/Tpetra specific code to access the underlying map data
820 
821  // check whether we have a Tpetra based Thyra operator
822  bool bIsTpetra = false;
823 #ifdef HAVE_XPETRA_TPETRA
824 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
825  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
826  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vectorSpace);
827  bIsTpetra = Teuchos::is_null(tpetra_vsc) ? false : true;
828 #endif
829 #endif
830 
831  // check whether we have an Epetra based Thyra operator
832  bool bIsEpetra = !bIsTpetra; // note: this is a little bit fragile!
833 
834 #ifdef HAVE_XPETRA_TPETRA
835  if(bIsTpetra) {
836 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
837  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
838  typedef Thyra::VectorBase<Scalar> ThyVecBase;
839  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
840  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
841  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
842  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
843  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
844  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
845  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
846  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
847  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
848 
849  resultMap = Xpetra::toXpetraNonConst(rgTpetraMap);
850 #else
851  throw Xpetra::Exceptions::RuntimeError("Problem AAA. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
852 #endif
853  }
854 #endif
855 
856 #ifdef HAVE_XPETRA_EPETRA
857  if(bIsEpetra) {
858  //RCP<const Epetra_Map> epMap = Teuchos::null;
859  RCP<const Epetra_Map>
860  epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map> >(vectorSpace,"epetra_map");
861  if(!Teuchos::is_null(epetra_map)) {
862  resultMap = Teuchos::rcp(new Xpetra::EpetraMapT<GlobalOrdinal,Node>(epetra_map));
863  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
864  } else {
865  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "No Epetra_Map data found in Thyra::VectorSpace.");
866  }
867  }
868 #endif
869  } // end standard case (no product map)
870  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
871  return resultMap;
872  }
873 
874  // const version
875  static Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
876  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
877  using Teuchos::RCP;
878  using Teuchos::rcp_dynamic_cast;
879  using Teuchos::as;
880 
881  using ThyProdMultVecBase = Thyra::ProductMultiVectorBase<Scalar>;
882  using ThyMultVecBase = Thyra::MultiVectorBase<Scalar>;
883  using ThyUtils = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
884 
885  // return value
886  RCP<MultiVector> xpMultVec = Teuchos::null;
887 
888  // check whether v is a product multi vector
889  Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase >(v);
890  if(thyProdVec != Teuchos::null) {
891  // SPECIAL CASE: create a nested BlockedMultiVector
892  // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
893  RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
894  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
895 
896  // create new Xpetra::BlockedMultiVector
897  xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
898 
899  RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec);
900  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpBlockedMultVec));
901 
902  // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
903  for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
904  RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
905  // xpBlockMV can be of type MultiVector or BlockedMultiVector
906  RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); //recursive call
907  xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
908  }
909 
910  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
911  return xpMultVec;
912  } else {
913  // STANDARD CASE: no product vector
914  // Epetra/Tpetra specific code to access the underlying map data
915  bool bIsTpetra = false;
916 #ifdef HAVE_XPETRA_TPETRA
917 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
918  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
919 
920  //typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
921  //typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
922  typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
923  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
924  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
925  typedef Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpTpMultVec;
926  typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
927 
928  RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
929  RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
930  if(thyraTpetraMultiVector != Teuchos::null) {
931  bIsTpetra = true;
932  const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
933  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
934  RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
935  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
936  xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
937  }
938 #endif
939 #endif
940 
941 #ifdef HAVE_XPETRA_EPETRA
942  if(bIsTpetra == false) {
943  // no product vector
944  Teuchos::RCP<Map> map = ThyUtils::toXpetra(v->range(), comm);
945  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(map));
946  RCP<Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = rcp_dynamic_cast<Xpetra::EpetraMapT<GlobalOrdinal,Node> >(map);
947  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xeMap));
948  RCP<const Epetra_Map> eMap = Teuchos::rcpFromRef(xeMap->getEpetra_Map());
949  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(eMap));
950  Teuchos::RCP<const Epetra_MultiVector> epMultVec = Thyra::get_Epetra_MultiVector(*eMap, v);
951  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epMultVec));
952  RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<Epetra_MultiVector>(epMultVec);
953  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epNonConstMultVec));
954  xpMultVec = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node>(epNonConstMultVec));
955  }
956 #endif
957  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
958  return xpMultVec;
959  } // end standard case
960  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
961  }
962 
963  // non-const version
964  static Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
965  toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
966  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > cv =
967  Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar> >(v);
968  Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > r =
969  toXpetra(cv,comm);
970  return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(r);
971  }
972 
973  static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
974  // check whether we have a Tpetra based Thyra operator
975  bool bIsTpetra = false;
976 #ifdef HAVE_XPETRA_TPETRA
977 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
978  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
979 
980  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(op);
981  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
982 
983  // for debugging purposes: find out why dynamic cast failed
984  if(!bIsTpetra &&
985 #ifdef HAVE_XPETRA_EPETRA
986  Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
987 #endif
988  Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
989  // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
990  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
991  if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
992  std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
993  std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
994  std::cout << " properly set!" << std::endl;
995  std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
996  }
997  }
998 #endif
999 #endif
1000 
1001 #if 0
1002  // Check whether it is a blocked operator.
1003  // If yes, grab the (0,0) block and check the underlying linear algebra
1004  // Note: we expect that the (0,0) block exists!
1005  if(bIsTpetra == false) {
1006  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1007  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
1008  if(ThyBlockedOp != Teuchos::null) {
1009  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1010  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1011  ThyBlockedOp->getBlock(0,0);
1012  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1013  bIsTpetra = isTpetra(b00);
1014  }
1015  }
1016 #endif
1017 
1018  return bIsTpetra;
1019  }
1020 
1021  static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
1022  // check whether we have an Epetra based Thyra operator
1023  bool bIsEpetra = false;
1024 
1025 #ifdef HAVE_XPETRA_EPETRA
1026  Teuchos::RCP<const Thyra::EpetraLinearOp> epetra_op = Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op,false);
1027  bIsEpetra = Teuchos::is_null(epetra_op) ? false : true;
1028 #endif
1029 
1030 #if 0
1031  // Check whether it is a blocked operator.
1032  // If yes, grab the (0,0) block and check the underlying linear algebra
1033  // Note: we expect that the (0,0) block exists!
1034  if(bIsEpetra == false) {
1035  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1036  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op,false);
1037  if(ThyBlockedOp != Teuchos::null) {
1038  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1039  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1040  ThyBlockedOp->getBlock(0,0);
1041  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1042  bIsEpetra = isEpetra(b00);
1043  }
1044  }
1045 #endif
1046 
1047  return bIsEpetra;
1048  }
1049 
1050  static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
1051  // Check whether it is a blocked operator.
1052  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1053  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
1054  if(ThyBlockedOp != Teuchos::null) {
1055  return true;
1056  }
1057  return false;
1058  }
1059 
1060  static Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
1061  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >& op) {
1062 
1063 #ifdef HAVE_XPETRA_TPETRA
1064  if(isTpetra(op)) {
1065 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1066  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1067 
1068  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1069  Teuchos::RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getConstTpetraOperator(op);
1070  // we should also add support for the const versions!
1071  //getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
1072  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1073  Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
1074  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraRowMat));
1075  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
1076  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraCrsMat));
1077  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraCrsMat);
1078  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
1079 
1080  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xTpetraCrsMat =
1081  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(TpetraNcnstCrsMat));
1082  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1083  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ret =
1084  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xTpetraCrsMat);
1085  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ret));
1086  return ret;
1087 #else
1088  throw Xpetra::Exceptions::RuntimeError("Problem BBB. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1089 #endif
1090  }
1091 #endif
1092 
1093 #ifdef HAVE_XPETRA_EPETRA
1094  if(isEpetra(op)) {
1095  Teuchos::RCP<const Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator( *op );
1096  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1097  Teuchos::RCP<const Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_op);
1098  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_rowmat));
1099  Teuchos::RCP<const Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(epetra_rowmat);
1100  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_crsmat));
1101  Teuchos::RCP<Epetra_CrsMatrix> epetra_ncnstcrsmat = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(epetra_crsmat);
1102  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_ncnstcrsmat));
1103 
1104  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > xEpetraCrsMat =
1105  Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> (epetra_ncnstcrsmat));
1106  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1107 
1108  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ret =
1109  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xEpetraCrsMat);
1110  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ret));
1111  return ret;
1112  }
1113 #endif
1114  return Teuchos::null;
1115  }
1116 
1117  static Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
1118  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar> >& op) {
1119 
1120 #ifdef HAVE_XPETRA_TPETRA
1121  if(isTpetra(op)) {
1122 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1123  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1124 
1125  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1126  Teuchos::RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getTpetraOperator(op);
1127  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1128  Teuchos::RCP<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
1129  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraRowMat));
1130  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
1131  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraCrsMat));
1132 
1133  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xTpetraCrsMat =
1134  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(TpetraCrsMat));
1135  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1136  return xTpetraCrsMat;
1137 #else
1138  throw Xpetra::Exceptions::RuntimeError("Problem CCC. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1139 #endif
1140  }
1141 #endif
1142 
1143 #ifdef HAVE_XPETRA_EPETRA
1144  if(isEpetra(op)) {
1145  Teuchos::RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator( *op );
1146  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1147  Teuchos::RCP<Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetra_op);
1148  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_rowmat));
1149  Teuchos::RCP<Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetra_rowmat);
1150  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_crsmat));
1151 
1152  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > xEpetraCrsMat =
1153  Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> (epetra_crsmat));
1154  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1155  return xEpetraCrsMat;
1156  }
1157 #endif
1158  return Teuchos::null;
1159  }
1160 
1161  static Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
1162  toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map) {
1163  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > thyraMap = Teuchos::null;
1164 
1165  // check whether map is of type BlockedMap
1166  RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
1167  if(bmap.is_null() == false) {
1168 
1169  Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > > vecSpaces(bmap->getNumMaps());
1170  for(size_t i = 0; i < bmap->getNumMaps(); i++) {
1171  // we need Thyra GIDs for all the submaps
1172  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > vs =
1173  Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(bmap->getMap(i,true));
1174  vecSpaces[i] = vs;
1175  }
1176 
1177  thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
1178  return thyraMap;
1179  }
1180 
1181  // standard case
1182 #ifdef HAVE_XPETRA_TPETRA
1183  if(map->lib() == Xpetra::UseTpetra) {
1184 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1185  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1186  Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> > tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(map);
1187  if (tpetraMap == Teuchos::null)
1188  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
1189  RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > tpMap = tpetraMap->getTpetra_Map();
1190  RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
1191  thyraMap = thyraTpetraMap;
1192 #else
1193  throw Xpetra::Exceptions::RuntimeError("Problem DDD. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1194 #endif
1195  }
1196 #endif
1197 
1198 #ifdef HAVE_XPETRA_EPETRA
1199  if(map->lib() == Xpetra::UseEpetra) {
1200  Teuchos::RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > epetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(map);
1201  if (epetraMap == Teuchos::null)
1202  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
1203  RCP<const Thyra::VectorSpaceBase<Scalar> > thyraEpetraMap = Thyra::create_VectorSpace(Teuchos::rcpFromRef(epetraMap->getEpetra_Map()));
1204  thyraMap = thyraEpetraMap;
1205  }
1206 #endif
1207 
1208  return thyraMap;
1209  }
1210 
1211  static Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
1212  toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > vec) {
1213 
1214  // create Thyra vector space out of Xpetra Map
1215  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > thMap = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(vec->getMap());
1216  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.");
1217  Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> > thSpmdMap = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<Scalar> >(thMap);
1218  TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error, "Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
1219  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.");
1220 
1221  // create Thyra MultiVector
1222  Teuchos::RCP< Thyra::MultiVectorBase<Scalar> > thMVec = Thyra::createMembers(thMap, vec->getNumVectors());
1223  Teuchos::RCP< Thyra::SpmdMultiVectorBase<Scalar> > thSpmdMVec = Teuchos::rcp_dynamic_cast<Thyra::SpmdMultiVectorBase<Scalar> >(thMVec);
1224  TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMVec == Teuchos::null, std::logic_error, "Cannot cast MultiVectorBase to SpmdMultiVectorBase.");
1225 
1226  // fill multivector with some data
1227  const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
1228  const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
1229  Teuchos::RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1230  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thSpmdMVec,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
1231 
1232  // loop over all vectors in multivector
1233  for(size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
1234  Teuchos::ArrayRCP< const Scalar > vecData = vec->getData(j);
1235  // loop over all local rows
1236  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1237  (*thyData)(i,j) = vecData[i];
1238  }
1239  }
1240 
1241  return thMVec;
1242  }
1243 
1244  static Teuchos::RCP<const Thyra::VectorBase<Scalar> >
1245  toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > vec) {
1246 
1247  // create Thyra vector space out of Xpetra Map
1248  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > thMap = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(vec->getMap());
1249  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.");
1250  Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> > thSpmdMap = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<Scalar> >(thMap);
1251  TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error, "Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
1252  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.");
1253 
1254  // create Thyra MultiVector
1255  Teuchos::RCP< Thyra::VectorBase<Scalar> > thMVec = Thyra::createMember(thMap);
1256  Teuchos::RCP< Thyra::SpmdVectorBase<Scalar> > thSpmdMVec = Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<Scalar> >(thMVec);
1257  TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMVec == Teuchos::null, std::logic_error, "Cannot cast VectorBase to SpmdVectorBase.");
1258 
1259  // fill multivector with some data
1260  const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
1261  const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
1262  Teuchos::RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1263  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thSpmdMVec,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
1264 
1265  // loop over all vectors in multivector
1266  for(size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
1267  Teuchos::ArrayRCP< const Scalar > vecData = vec->getData(j);
1268  // loop over all local rows
1269  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1270  (*thyData)(i,j) = vecData[i];
1271  }
1272  }
1273 
1274  return thMVec;
1275  }
1276 
1277  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) {
1278  using Teuchos::RCP;
1279  using Teuchos::rcp_dynamic_cast;
1280  using Teuchos::as;
1281  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1282  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
1283  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1284  //typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
1285  //typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
1286  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1287 
1288  // copy data from tY_inout to Y_inout
1289  RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
1290  if(prodTarget != Teuchos::null) {
1291 
1292  RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
1293  if(bSourceVec.is_null() == true) {
1294  // SPECIAL CASE: target vector is product vector:
1295  // update Thyra product multi vector with data from a merged Xpetra multi vector
1296 
1297  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
1298  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.");
1299 
1300  for(int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1301  // access Xpetra data
1302  RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
1303 
1304  // access Thyra data
1305  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
1306  RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
1307  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
1308  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
1309  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
1310  RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1311  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
1312 
1313  // loop over all vectors in multivector
1314  for(size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
1315  Teuchos::ArrayRCP< const Scalar > xpData = xpSubBlock->getData(j); // access const data from Xpetra object
1316 
1317  // loop over all local rows
1318  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1319  (*thyData)(i,j) = xpData[i];
1320  }
1321  }
1322  }
1323  } else {
1324  // source vector is a blocked multivector
1325  // TODO test me
1326  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.");
1327 
1328  for(int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1329  // access Thyra data
1330  RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
1331 
1332  Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
1333 
1334  // access Thyra data
1335  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
1336  Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
1337  }
1338 
1339  }
1340  } else {
1341  // STANDARD case:
1342  // update Thyra::MultiVector with data from an Xpetra::MultiVector
1343 
1344  // access Thyra data
1345  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
1346  TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
1347  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
1348  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
1349  RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1350  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
1351 
1352  // loop over all vectors in multivector
1353  for(size_t j = 0; j < source->getNumVectors(); ++j) {
1354  Teuchos::ArrayRCP< const Scalar > xpData = source->getData(j); // access const data from Xpetra object
1355  // loop over all local rows
1356  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1357  (*thyData)(i,j) = xpData[i];
1358  }
1359  }
1360  }
1361  }
1362 
1363  static Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
1364  toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& mat) {
1365  // create a Thyra operator from Xpetra::CrsMatrix
1366  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
1367 
1368 #ifdef HAVE_XPETRA_TPETRA
1369  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
1370  if(tpetraMat!=Teuchos::null) {
1371 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1372  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1373 
1374  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
1375  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpCrsMat));
1376  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
1377  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
1378 
1379  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
1380  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpRowMat));
1381  Teuchos::RCP<const Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
1382  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpOperator));
1383 
1384  thyraOp = Thyra::createConstLinearOp(tpOperator);
1385  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
1386 #else
1387  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1388 #endif
1389  }
1390 #endif
1391 
1392 #ifdef HAVE_XPETRA_EPETRA
1393  Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > epetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(mat);
1394  if(epetraMat!=Teuchos::null) {
1395  Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > xEpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(mat);
1396  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpCrsMat));
1397  Teuchos::RCP<const Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrix();
1398  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
1399 
1400  Teuchos::RCP<const Thyra::EpetraLinearOp> thyraEpOp = Thyra::epetraLinearOp(epCrsMat,"op");
1401  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
1402  thyraOp = thyraEpOp;
1403  }
1404 #endif
1405  return thyraOp;
1406  }
1407 
1408  static Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1409  toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& mat) {
1410  // create a Thyra operator from Xpetra::CrsMatrix
1411  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
1412 
1413 #ifdef HAVE_XPETRA_TPETRA
1414  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
1415  if(tpetraMat!=Teuchos::null) {
1416 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1417  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1418 
1419  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
1420  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpCrsMat));
1421  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
1422  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
1423 
1424  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
1425  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpRowMat));
1426  Teuchos::RCP<Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
1427  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpOperator));
1428 
1429  thyraOp = Thyra::createLinearOp(tpOperator);
1430  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
1431 #else
1432  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1433 #endif
1434  }
1435 #endif
1436 
1437 #ifdef HAVE_XPETRA_EPETRA
1438  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(mat);
1439  if(epetraMat!=Teuchos::null) {
1440  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > xEpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(mat);
1441  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpCrsMat));
1442  Teuchos::RCP<Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrixNonConst();
1443  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
1444 
1445  Teuchos::RCP<Thyra::EpetraLinearOp> thyraEpOp = Thyra::nonconstEpetraLinearOp(epCrsMat,"op");
1446  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
1447  thyraOp = thyraEpOp;
1448  }
1449 #endif
1450  return thyraOp;
1451  }
1452 
1453  static Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
1454  toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, int, EpetraNode> >& mat);
1455 
1456 }; // specialization on SC=double, LO=GO=int and NO=EpetraNode
1457 #endif // #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
1458 
1459 #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
1460 template <>
1461 class ThyraUtils<double, int, long long, EpetraNode> {
1462 public:
1463  typedef double Scalar;
1464  typedef int LocalOrdinal;
1465  typedef long long GlobalOrdinal;
1466  typedef EpetraNode Node;
1467 
1468 private:
1469 #undef XPETRA_THYRAUTILS_SHORT
1470 #include "Xpetra_UseShortNames.hpp"
1471 
1472 public:
1473 
1474  static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >
1475  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) {
1476 
1477  Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map = ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(vectorSpace,comm);
1478 
1479  if(stridedBlockId == -1) {
1480  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo.size() != 0);
1481  }
1482  else {
1483  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo[stridedBlockId] != 0);
1484  }
1485 
1486  Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> > ret = Xpetra::StridedMapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(map, stridingInfo, stridedBlockId, offset);
1487  return ret;
1488  }
1489 
1490  static Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
1491  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
1492  using Teuchos::RCP;
1493  using Teuchos::rcp_dynamic_cast;
1494  using Teuchos::as;
1495  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1496  typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
1497  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
1498 
1499  RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase >(vectorSpace);
1500  if(prodVectorSpace != Teuchos::null) {
1501  // SPECIAL CASE: product Vector space
1502  // collect all submaps to store them in a hierarchical BlockedMap object
1503  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error, "Found a product vector space with zero blocks.");
1504  std::vector<RCP<const Map> > mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
1505  std::vector<RCP<const Map> > mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
1506  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
1507  RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
1508  // can be of type Map or BlockedMap (containing Thyra GIDs)
1509  mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
1510  }
1511 
1512  // get offsets for submap GIDs
1513  // we need that for the full map (Xpetra GIDs)
1514  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
1515  for(int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
1516  gidOffsets[i] = mapsThyra[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
1517  }
1518 
1519  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
1520  RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
1521  // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
1522  mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
1523  }
1524 
1525  Teuchos::RCP<Map> resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(mapsXpetra, mapsThyra));
1526  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
1527  return resultMap;
1528  } else {
1529  // STANDARD CASE: no product map
1530  // Epetra/Tpetra specific code to access the underlying map data
1531 
1532  // check whether we have a Tpetra based Thyra operator
1533  bool bIsTpetra = false;
1534 #ifdef HAVE_XPETRA_TPETRA
1535 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1536  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1537  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vectorSpace);
1538  bIsTpetra = Teuchos::is_null(tpetra_vsc) ? false : true;
1539 #endif
1540 #endif
1541 
1542  // check whether we have an Epetra based Thyra operator
1543  bool bIsEpetra = !bIsTpetra; // note: this is a little bit fragile!
1544 
1545 #ifdef HAVE_XPETRA_TPETRA
1546  if(bIsTpetra) {
1547 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1548  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1549  typedef Thyra::VectorBase<Scalar> ThyVecBase;
1550  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
1551  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
1552  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1553  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
1554  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
1555  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
1556  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
1557  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
1558  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
1559 
1560  RCP<Map> rgXpetraMap = Xpetra::toXpetraNonConst(rgTpetraMap);
1561  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgXpetraMap));
1562  return rgXpetraMap;
1563 #else
1564  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1565 #endif
1566  }
1567 #endif
1568 
1569 #ifdef HAVE_XPETRA_EPETRA
1570  if(bIsEpetra) {
1571  //RCP<const Epetra_Map> epMap = Teuchos::null;
1572  RCP<const Epetra_Map>
1573  epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map> >(vectorSpace,"epetra_map");
1574  if(!Teuchos::is_null(epetra_map)) {
1575  Teuchos::RCP<Map> rgXpetraMap = Teuchos::rcp(new Xpetra::EpetraMapT<GlobalOrdinal,Node>(epetra_map));
1576  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgXpetraMap));
1577  return rgXpetraMap;
1578  } else {
1579  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "No Epetra_Map data found in Thyra::VectorSpace.");
1580  }
1581  }
1582 #endif
1583  } // end standard case (no product map)
1584  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "Cannot transform Thyra::VectorSpace to Xpetra::Map.");
1585  // return Teuchos::null; // unreachable
1586  }
1587 
1588  // const version
1589  static Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
1590  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
1591  using Teuchos::RCP;
1592  using Teuchos::rcp_dynamic_cast;
1593  using Teuchos::as;
1594  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1595  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1596  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
1597 
1598  // return value
1599  RCP<MultiVector> xpMultVec = Teuchos::null;
1600 
1601  // check whether v is a product multi vector
1602  Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase >(v);
1603  if(thyProdVec != Teuchos::null) {
1604  // SPECIAL CASE: create a nested BlockedMultiVector
1605  // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
1606  RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
1607  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
1608 
1609  // create new Xpetra::BlockedMultiVector
1610  xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
1611 
1612  RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec);
1613  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpBlockedMultVec));
1614 
1615  // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
1616  for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
1617  RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
1618  // xpBlockMV can be of type MultiVector or BlockedMultiVector
1619  RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); //recursive call
1620  xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
1621  }
1622 
1623  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1624  return xpMultVec;
1625  } else {
1626  // STANDARD CASE: no product vector
1627  // Epetra/Tpetra specific code to access the underlying map data
1628  bool bIsTpetra = false;
1629 #ifdef HAVE_XPETRA_TPETRA
1630 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1631  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1632 
1633  //typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
1634  //typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
1635  typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
1636  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
1637  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
1638  typedef Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpTpMultVec;
1639  typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
1640 
1641  RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
1642  RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
1643  if(thyraTpetraMultiVector != Teuchos::null) {
1644  bIsTpetra = true;
1645  const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
1646  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
1647  RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
1648  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
1649  xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
1650  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1651  return xpMultVec;
1652  }
1653 #endif
1654 #endif
1655 
1656 #ifdef HAVE_XPETRA_EPETRA
1657  if(bIsTpetra == false) {
1658  // no product vector
1659  Teuchos::RCP<Map> map = ThyUtils::toXpetra(v->range(), comm);
1660  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(map));
1661  RCP<Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = rcp_dynamic_cast<Xpetra::EpetraMapT<GlobalOrdinal,Node> >(map);
1662  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xeMap));
1663  RCP<const Epetra_Map> eMap = Teuchos::rcpFromRef(xeMap->getEpetra_Map());
1664  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(eMap));
1665  Teuchos::RCP<const Epetra_MultiVector> epMultVec = Thyra::get_Epetra_MultiVector(*eMap, v);
1666  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epMultVec));
1667  RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<Epetra_MultiVector>(epMultVec);
1668  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epNonConstMultVec));
1669  xpMultVec = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node>(epNonConstMultVec));
1670  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
1671  return xpMultVec;
1672  }
1673 #endif
1674  } // end standard case
1675  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "Cannot transform Thyra::MultiVector to Xpetra::MultiVector.");
1676  // return Teuchos::null; // unreachable
1677  }
1678 
1679  // non-const version
1680  static Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
1681  toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
1682  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > cv =
1683  Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar> >(v);
1684  Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > r =
1685  toXpetra(cv,comm);
1686  return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(r);
1687  }
1688 
1689  static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
1690  // check whether we have a Tpetra based Thyra operator
1691  bool bIsTpetra = false;
1692 #ifdef HAVE_XPETRA_TPETRA
1693 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1694  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1695 
1696  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(op);
1697  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
1698 
1699  // for debugging purposes: find out why dynamic cast failed
1700  if(!bIsTpetra &&
1701 #ifdef HAVE_XPETRA_EPETRA
1702  Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
1703 #endif
1704  Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
1705  // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
1706  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
1707  if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
1708  std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
1709  std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
1710  std::cout << " properly set!" << std::endl;
1711  std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
1712  }
1713  }
1714 #endif
1715 #endif
1716 
1717 #if 0
1718  // Check whether it is a blocked operator.
1719  // If yes, grab the (0,0) block and check the underlying linear algebra
1720  // Note: we expect that the (0,0) block exists!
1721  if(bIsTpetra == false) {
1722  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1723  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
1724  if(ThyBlockedOp != Teuchos::null) {
1725  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1726  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1727  ThyBlockedOp->getBlock(0,0);
1728  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1729  bIsTpetra = isTpetra(b00);
1730  }
1731  }
1732 #endif
1733 
1734  return bIsTpetra;
1735  }
1736 
1737  static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
1738  // check whether we have an Epetra based Thyra operator
1739  bool bIsEpetra = false;
1740 
1741 #ifdef HAVE_XPETRA_EPETRA
1742  Teuchos::RCP<const Thyra::EpetraLinearOp> epetra_op = Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op,false);
1743  bIsEpetra = Teuchos::is_null(epetra_op) ? false : true;
1744 #endif
1745 
1746 #if 0
1747  // Check whether it is a blocked operator.
1748  // If yes, grab the (0,0) block and check the underlying linear algebra
1749  // Note: we expect that the (0,0) block exists!
1750  if(bIsEpetra == false) {
1751  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1752  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op,false);
1753  if(ThyBlockedOp != Teuchos::null) {
1754  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
1755  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
1756  ThyBlockedOp->getBlock(0,0);
1757  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
1758  bIsEpetra = isEpetra(b00);
1759  }
1760  }
1761 #endif
1762 
1763  return bIsEpetra;
1764  }
1765 
1766  static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
1767  // Check whether it is a blocked operator.
1768  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
1769  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
1770  if(ThyBlockedOp != Teuchos::null) {
1771  return true;
1772  }
1773  return false;
1774  }
1775 
1776  static Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
1777  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >& op) {
1778 
1779 #ifdef HAVE_XPETRA_TPETRA
1780  if(isTpetra(op)) {
1781 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1782  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1783 
1784  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1785  Teuchos::RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getConstTpetraOperator(op);
1786  // we should also add support for the const versions!
1787  //getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
1788  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1789  Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
1790  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraRowMat));
1791  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
1792  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraCrsMat));
1793  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraCrsMat);
1794  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
1795 
1796  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xTpetraCrsMat =
1797  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(TpetraNcnstCrsMat));
1798  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1799  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ret =
1800  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xTpetraCrsMat);
1801  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ret));
1802  return ret;
1803 #else
1804  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1805 #endif
1806  }
1807 #endif
1808 
1809 #ifdef HAVE_XPETRA_EPETRA
1810  if(isEpetra(op)) {
1811  Teuchos::RCP<const Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator( *op );
1812  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1813  Teuchos::RCP<const Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_op);
1814  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_rowmat));
1815  Teuchos::RCP<const Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(epetra_rowmat);
1816  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_crsmat));
1817  Teuchos::RCP<Epetra_CrsMatrix> epetra_ncnstcrsmat = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(epetra_crsmat);
1818  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_ncnstcrsmat));
1819 
1820  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > xEpetraCrsMat =
1821  Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> (epetra_ncnstcrsmat));
1822  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1823 
1824  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ret =
1825  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xEpetraCrsMat);
1826  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ret));
1827  return ret;
1828  }
1829 #endif
1830  return Teuchos::null;
1831  }
1832 
1833  static Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
1834  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar> >& op) {
1835 
1836 #ifdef HAVE_XPETRA_TPETRA
1837  if(isTpetra(op)) {
1838 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1839  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1840 
1841  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1842  Teuchos::RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getTpetraOperator(op);
1843  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
1844  Teuchos::RCP<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
1845  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraRowMat));
1846  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
1847  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraCrsMat));
1848 
1849  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xTpetraCrsMat =
1850  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(TpetraCrsMat));
1851  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1852  return xTpetraCrsMat;
1853 #else
1854  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1855 #endif
1856  }
1857 #endif
1858 
1859 #ifdef HAVE_XPETRA_EPETRA
1860  if(isEpetra(op)) {
1861  Teuchos::RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator( *op );
1862  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
1863  Teuchos::RCP<Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetra_op);
1864  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_rowmat));
1865  Teuchos::RCP<Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetra_rowmat);
1866  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_crsmat));
1867 
1868  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > xEpetraCrsMat =
1869  Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> (epetra_crsmat));
1870  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1871  return xEpetraCrsMat;
1872  }
1873 #endif
1874  return Teuchos::null;
1875  }
1876 
1877  static Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
1878  toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map) {
1879  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > thyraMap = Teuchos::null;
1880 
1881  // check whether map is of type BlockedMap
1882  RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
1883  if(bmap.is_null() == false) {
1884 
1885  Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > > vecSpaces(bmap->getNumMaps());
1886  for(size_t i = 0; i < bmap->getNumMaps(); i++) {
1887  // we need Thyra GIDs for all the submaps
1888  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > vs =
1889  Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(bmap->getMap(i,true));
1890  vecSpaces[i] = vs;
1891  }
1892 
1893  thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
1894  return thyraMap;
1895  }
1896 
1897  // standard case
1898 #ifdef HAVE_XPETRA_TPETRA
1899  if(map->lib() == Xpetra::UseTpetra) {
1900 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1901  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1902  Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> > tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(map);
1903  if (tpetraMap == Teuchos::null)
1904  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
1905  RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > tpMap = tpetraMap->getTpetra_Map();
1906  RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
1907  thyraMap = thyraTpetraMap;
1908 #else
1909  throw Xpetra::Exceptions::RuntimeError("Problem DDD. Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
1910 #endif
1911  }
1912 #endif
1913 
1914 #ifdef HAVE_XPETRA_EPETRA
1915  if(map->lib() == Xpetra::UseEpetra) {
1916  Teuchos::RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > epetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(map);
1917  if (epetraMap == Teuchos::null)
1918  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
1919  RCP<const Thyra::VectorSpaceBase<Scalar> > thyraEpetraMap = Thyra::create_VectorSpace(Teuchos::rcpFromRef(epetraMap->getEpetra_Map()));
1920  thyraMap = thyraEpetraMap;
1921  }
1922 #endif
1923 
1924  return thyraMap;
1925  }
1926 
1927  static Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
1928  toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > vec) {
1929 
1930  // create Thyra vector space out of Xpetra Map
1931  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > thMap = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(vec->getMap());
1932  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.");
1933  Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> > thSpmdMap = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<Scalar> >(thMap);
1934  TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error, "Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
1935  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.");
1936 
1937  // create Thyra MultiVector
1938  Teuchos::RCP< Thyra::MultiVectorBase<Scalar> > thMVec = Thyra::createMembers(thMap, vec->getNumVectors());
1939  Teuchos::RCP< Thyra::SpmdMultiVectorBase<Scalar> > thSpmdMVec = Teuchos::rcp_dynamic_cast<Thyra::SpmdMultiVectorBase<Scalar> >(thMVec);
1940  TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMVec == Teuchos::null, std::logic_error, "Cannot cast MultiVectorBase to SpmdMultiVectorBase.");
1941 
1942  // fill multivector with some data
1943  const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
1944  const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
1945  Teuchos::RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1946  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thSpmdMVec,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
1947 
1948  // loop over all vectors in multivector
1949  for(size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
1950  Teuchos::ArrayRCP< const Scalar > vecData = vec->getData(j);
1951  // loop over all local rows
1952  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1953  (*thyData)(i,j) = vecData[i];
1954  }
1955  }
1956 
1957  return thMVec;
1958  }
1959 
1960  static Teuchos::RCP<const Thyra::VectorBase<Scalar> >
1961  toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > vec) {
1962 
1963  // create Thyra vector space out of Xpetra Map
1964  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > thMap = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(vec->getMap());
1965  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.");
1966  Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<Scalar> > thSpmdMap = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<Scalar> >(thMap);
1967  TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error, "Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
1968  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.");
1969 
1970  // create Thyra MultiVector
1971  Teuchos::RCP< Thyra::VectorBase<Scalar> > thMVec = Thyra::createMember(thMap);
1972  Teuchos::RCP< Thyra::SpmdVectorBase<Scalar> > thSpmdMVec = Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<Scalar> >(thMVec);
1973  TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMVec == Teuchos::null, std::logic_error, "Cannot cast VectorBase to SpmdVectorBase.");
1974 
1975  // fill multivector with some data
1976  const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
1977  const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
1978  Teuchos::RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1979  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thSpmdMVec,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
1980 
1981  // loop over all vectors in multivector
1982  for(size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
1983  Teuchos::ArrayRCP< const Scalar > vecData = vec->getData(j);
1984  // loop over all local rows
1985  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1986  (*thyData)(i,j) = vecData[i];
1987  }
1988  }
1989 
1990  return thMVec;
1991  }
1992 
1993  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) {
1994  using Teuchos::RCP;
1995  using Teuchos::rcp_dynamic_cast;
1996  using Teuchos::as;
1997  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1998  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
1999  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
2000  //typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
2001  //typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
2002  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
2003 
2004  // copy data from tY_inout to Y_inout
2005  RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
2006  if(prodTarget != Teuchos::null) {
2007  RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
2008  if(bSourceVec.is_null() == true) {
2009  // SPECIAL CASE: target vector is product vector:
2010  // update Thyra product multi vector with data from a merged Xpetra multi vector
2011 
2012  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
2013  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.");
2014 
2015  for(int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
2016  // access Xpetra data
2017  RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
2018 
2019  // access Thyra data
2020  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
2021  RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
2022  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
2023  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
2024  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
2025  RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
2026  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
2027 
2028  // loop over all vectors in multivector
2029  for(size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
2030  Teuchos::ArrayRCP< const Scalar > xpData = xpSubBlock->getData(j); // access const data from Xpetra object
2031 
2032  // loop over all local rows
2033  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
2034  (*thyData)(i,j) = xpData[i];
2035  }
2036  }
2037  }
2038  } else {
2039  // source vector is a blocked multivector
2040  // TODO test me
2041  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.");
2042 
2043  for(int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
2044  // access Thyra data
2045  RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
2046 
2047  Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
2048 
2049  // access Thyra data
2050  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
2051  Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
2052  }
2053 
2054  }
2055  } else {
2056  // STANDARD case:
2057  // update Thyra::MultiVector with data from an Xpetra::MultiVector
2058 
2059  // access Thyra data
2060  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
2061  TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
2062  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
2063  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
2064  RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
2065  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
2066 
2067  // loop over all vectors in multivector
2068  for(size_t j = 0; j < source->getNumVectors(); ++j) {
2069  Teuchos::ArrayRCP< const Scalar > xpData = source->getData(j); // access const data from Xpetra object
2070  // loop over all local rows
2071  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
2072  (*thyData)(i,j) = xpData[i];
2073  }
2074  }
2075  }
2076  }
2077 
2078  static Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
2079  toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& mat) {
2080  // create a Thyra operator from Xpetra::CrsMatrix
2081  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
2082 
2083 #ifdef HAVE_XPETRA_TPETRA
2084  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
2085  if(tpetraMat!=Teuchos::null) {
2086 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
2087  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
2088 
2089  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
2090  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpCrsMat));
2091  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
2092  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
2093 
2094  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
2095  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpRowMat));
2096  Teuchos::RCP<const Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
2097  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpOperator));
2098 
2099  thyraOp = Thyra::createConstLinearOp(tpOperator);
2100  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
2101 #else
2102  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
2103 #endif
2104  }
2105 #endif
2106 
2107 #ifdef HAVE_XPETRA_EPETRA
2108  Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > epetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(mat);
2109  if(epetraMat!=Teuchos::null) {
2110  Teuchos::RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > xEpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(mat);
2111  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpCrsMat));
2112  Teuchos::RCP<const Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrix();
2113  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
2114 
2115  Teuchos::RCP<const Thyra::EpetraLinearOp> thyraEpOp = Thyra::epetraLinearOp(epCrsMat,"op");
2116  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
2117  thyraOp = thyraEpOp;
2118  }
2119 #endif
2120  return thyraOp;
2121  }
2122 
2123  static Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
2124  toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& mat) {
2125  // create a Thyra operator from Xpetra::CrsMatrix
2126  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
2127 
2128 #ifdef HAVE_XPETRA_TPETRA
2129  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
2130  if(tpetraMat!=Teuchos::null) {
2131 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
2132  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
2133 
2134  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
2135  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpCrsMat));
2136  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
2137  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
2138 
2139  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
2140  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpRowMat));
2141  Teuchos::RCP<Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
2142  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpOperator));
2143 
2144  thyraOp = Thyra::createLinearOp(tpOperator);
2145  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
2146 #else
2147  throw Xpetra::Exceptions::RuntimeError("Add TPETRA_INST_INT_LONG_LONG:BOOL=ON in your configuration.");
2148 #endif
2149  }
2150 #endif
2151 
2152 #ifdef HAVE_XPETRA_EPETRA
2153  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(mat);
2154  if(epetraMat!=Teuchos::null) {
2155  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > xEpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(mat);
2156  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpCrsMat));
2157  Teuchos::RCP<Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrixNonConst();
2158  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
2159 
2160  Teuchos::RCP<Thyra::EpetraLinearOp> thyraEpOp = Thyra::nonconstEpetraLinearOp(epCrsMat,"op");
2161  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
2162  thyraOp = thyraEpOp;
2163  }
2164 #endif
2165  return thyraOp;
2166  }
2167 
2168  static Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
2169  toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, long long, EpetraNode> >& mat);
2170 
2171 }; // specialization on SC=double, LO=GO=int and NO=EpetraNode
2172 #endif // XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
2173 
2174 #endif // HAVE_XPETRA_EPETRA
2175 
2176 } // end namespace Xpetra
2177 
2178 #define XPETRA_THYRAUTILS_SHORT
2179 #endif // HAVE_XPETRA_THYRA
2180 
2181 #endif // XPETRA_THYRAUTILS_HPP
Xpetra_ConfigDefs.hpp
Xpetra_EpetraCrsMatrix.hpp
Xpetra_EpetraMap.hpp
Xpetra_MapUtils.hpp
Xpetra_Matrix.hpp
Xpetra_UseShortNames.hpp
Teuchos::RCP::is_null
bool is_null() const
Teuchos::RCP::ptr
Ptr< T > ptr() const
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:102
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:232
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
Definition: Xpetra_MultiVectorFactory_def.hpp:62
Xpetra::StridedMapFactory::Build
static RCP< Xpetra::StridedMap< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
Definition: Xpetra_StridedMapFactory_def.hpp:77
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
Teuchos::as
TypeTo as(const TypeFrom &t)
TEUCHOS_UNREACHABLE_RETURN
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra
Xpetra namespace
Definition: Xpetra_BlockedCrsMatrix.hpp:88
Xpetra::toXpetraNonConst
const RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > toXpetraNonConst(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &map)
Definition: Xpetra_TpetraMap_decl.hpp:276
Xpetra::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Definition: Xpetra_EpetraCrsGraph.cpp:65