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