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