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