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