Xpetra_ThyraUtils.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef XPETRA_THYRAUTILS_HPP
48 #define XPETRA_THYRAUTILS_HPP
49 
50 #include "Xpetra_ConfigDefs.hpp"
51 #ifdef HAVE_XPETRA_THYRA
52 
53 #include <typeinfo>
54 
55 #ifdef HAVE_XPETRA_TPETRA
56 #include "Tpetra_ConfigDefs.hpp"
57 #endif
58 
59 #ifdef HAVE_XPETRA_EPETRA
60 #include "Epetra_config.h"
61 #include "Epetra_CombineMode.h"
62 #endif
63 
64 #include "Xpetra_Map.hpp"
65 #include "Xpetra_StridedMap.hpp"
67 #include "Xpetra_MapExtractor.hpp"
68 
69 #include <Thyra_VectorSpaceBase.hpp>
70 #include <Thyra_ProductVectorSpaceBase.hpp>
71 #include <Thyra_ProductMultiVectorBase.hpp>
72 #include <Thyra_VectorSpaceBase.hpp>
73 #include <Thyra_DefaultBlockedLinearOp.hpp>
74 #include <Thyra_LinearOpBase.hpp>
75 #include <Thyra_DetachedMultiVectorView.hpp>
76 
77 #ifdef HAVE_XPETRA_TPETRA
78 #include <Thyra_TpetraThyraWrappers.hpp>
79 #include <Thyra_TpetraVector.hpp>
80 #include <Thyra_TpetraVectorSpace.hpp>
81 #include <Tpetra_Map.hpp>
82 #include <Tpetra_Vector.hpp>
83 #include <Tpetra_CrsMatrix.hpp>
84 #include <Xpetra_TpetraMap.hpp>
86 #endif
87 #ifdef HAVE_XPETRA_EPETRA
88 #include <Thyra_EpetraLinearOp.hpp>
89 #include <Thyra_EpetraThyraWrappers.hpp>
90 #include <Thyra_SpmdVectorBase.hpp>
91 #include <Thyra_get_Epetra_Operator.hpp>
92 #include <Epetra_Map.h>
93 #include <Epetra_Vector.h>
94 #include <Epetra_CrsMatrix.h>
95 #include <Xpetra_EpetraMap.hpp>
97 #endif
98 
99 namespace Xpetra {
100 
101 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> class BlockedCrsMatrix;
102 
103 template <class Scalar,
104 class LocalOrdinal = int,
105 class GlobalOrdinal = LocalOrdinal,
107 class ThyraUtils {
108 
109 private:
110 #undef XPETRA_THYRAUTILS_SHORT
111 #include "Xpetra_UseShortNames.hpp"
112 
113 public:
114 
116  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) {
117 
119 
120  if(stridedBlockId == -1) {
121  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo.size() != 0);
122  }
123  else {
124  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo[stridedBlockId] != 0);
125  }
126 
128  return ret;
129  }
130 
132  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
133  using Teuchos::RCP;
134  using Teuchos::rcp_dynamic_cast;
135  using Teuchos::as;
136  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
137 //#ifdef HAVE_XPETRA_TPETRA
138 //#ifdef HAVE_XPETRA_TPETRA_INST_INT_INT
139 
140 //#endif
141 //#endif
142  typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
145  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
146 
147  RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase >(vectorSpace);
148  if(prodVectorSpace != Teuchos::null) {
149  // SPECIAL CASE: product Vector space
150  // merge all submaps to one large Xpetra Map
151  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error, "Found a product vector space with zero blocks.");
152  std::vector<RCP<Map> > maps(prodVectorSpace->numBlocks(), Teuchos::null);
153  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
154  RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
155  RCP<Map> map = ThyUtils::toXpetra(bv, comm); // recursive call
156  maps[b] = map;
157  }
158 
159  // get offsets for submap GIDs
160  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
161  for(int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
162  gidOffsets[i] = maps[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
163  }
164 
165  // loop over all sub maps and collect GIDs
166  std::vector<GlobalOrdinal> gids;
167  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
168  for(LocalOrdinal l = 0; l < as<LocalOrdinal>(maps[b]->getNodeNumElements()); ++l) {
169  GlobalOrdinal gid = maps[b]->getGlobalElement(l) + gidOffsets[b];
170  gids.push_back(gid);
171  }
172  }
173 
174  // create full map
176  Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
177  RCP<Map> fullMap = MapFactory::Build(maps[0]->lib(), INVALID, gidsView, maps[0]->getIndexBase(), comm);
178 
180  return fullMap;
181  } else {
182 #ifdef HAVE_XPETRA_TPETRA
183  // STANDARD CASE: no product map
184  // check whether we have a Tpetra based Thyra operator
185  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vectorSpace);
186  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.");
187 
188  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
189  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
190  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
191  typedef Thyra::VectorBase<Scalar> ThyVecBase;
192  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
194  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
196  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
198 
199  RCP<Map> rgXpetraMap = Xpetra::toXpetraNonConst(rgTpetraMap);
201  return rgXpetraMap;
202 #else
203  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.");
204 #endif
205  } // end standard case (no product map)
206  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "Cannot transform Thyra::VectorSpace to Xpetra::Map.");
207  return Teuchos::null;
208  }
209 
210  // const version
212  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
213  using Teuchos::RCP;
214  using Teuchos::rcp_dynamic_cast;
215  using Teuchos::as;
216  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
217  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
218  //typedef Thyra::VectorBase<Scalar> ThyVecBase;
219 //#ifdef HAVE_XPETRA_TPETRA
220 //#ifdef HAVE_XPETRA_TPETRA_INST_INT_INT
221 
222 //#endif
223 //#endif
224  //typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
225  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
227  //typedef Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node> MapFactory;
230  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
231  //typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
232  //typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
233 
234  // return value
235  RCP<MultiVector> xpMultVec = Teuchos::null;
236 
237  // check whether v is a product multi vector
238  Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase >(v);
239  if(thyProdVec != Teuchos::null) {
240  // SPECIAL CASE: product Vector
241  // merge all subvectors to one large Xpetra vector
242  RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
243 
244  // create new Epetra_MultiVector living on full map (stored in eMap)
245  xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
246 
247  // fill xpMultVec with Thyra data
248  std::vector<GlobalOrdinal> lidOffsets(thyProdVec->productSpace()->numBlocks()+1,0);
249  for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
250  Teuchos::RCP<const ThyVecSpaceBase> thySubMap = thyProdVec->productSpace()->getBlock(b);
251  Teuchos::RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(thySubMap);
252  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
253  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : thySubMap->dim() );
254  lidOffsets[b+1] = localSubDim + lidOffsets[b]; // calculate lid offset for next block
255  RCP<Thyra::ConstDetachedMultiVectorView<Scalar> > thyData =
256  Teuchos::rcp(new Thyra::ConstDetachedMultiVectorView<Scalar>(thyProdVec->getMultiVectorBlock(b),Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
257  for(size_t vv = 0; vv < xpMultVec->getNumVectors(); ++vv) {
258  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
259  xpMultVec->replaceLocalValue(i + lidOffsets[b] , vv, (*thyData)(i,vv));
260  }
261  }
262  }
264  return xpMultVec;
265  } else {
266  // STANDARD CASE: no product vector
267 #ifdef HAVE_XPETRA_TPETRA
268  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
269  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
271  typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
272  typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
273 
274  RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
275  RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
276  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.");
277  const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
279  RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
280  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
281  xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
283  return xpMultVec;
284 #else
285  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.");
286 #endif
287  } // end standard case
288  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "Cannot transform Thyra::MultiVector to Xpetra::MultiVector.");
289  return Teuchos::null;
290  }
291 
292  // non-const version
294  toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
296  Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar> >(v);
298  toXpetra(cv,comm);
299  return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(r);
300  }
301 
302 
303  static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
305 
306  // check whether we have a Tpetra based Thyra operator
307  bool bIsTpetra = false;
308 #ifdef HAVE_XPETRA_TPETRA
309  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(op);
310  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
311 
312  // for debugging purposes: find out why dynamic cast failed
313  if(!bIsTpetra &&
314 #ifdef HAVE_XPETRA_EPETRA
315  Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
316 #endif
317  Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
318  // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
319  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
320  if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
321  std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
322  std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
323  std::cout << " properly set!" << std::endl;
324  std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
325  }
326  }
327 #endif
328 
329 #if 0
330  // Check whether it is a blocked operator.
331  // If yes, grab the (0,0) block and check the underlying linear algebra
332  // Note: we expect that the (0,0) block exists!
333  if(bIsTpetra == false) {
335  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
336  if(ThyBlockedOp != Teuchos::null) {
337  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
339  ThyBlockedOp->getBlock(0,0);
341  bIsTpetra = isTpetra(b00);
342  }
343  }
344 #endif
345 
346  return bIsTpetra;
347  }
348 
349  static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
350  return false;
351  }
352 
353  static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
354  // Check whether it is a blocked operator.
356  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
357  if(ThyBlockedOp != Teuchos::null) {
358  return true;
359  }
360  return false;
361  }
362 
364  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >& op) {
365 
366 #ifdef HAVE_XPETRA_TPETRA
367  if(isTpetra(op)) {
368  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
369  Teuchos::RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getConstTpetraOperator(op);
370  // we should also add support for the const versions!
371  //getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
373  Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
375  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
377  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraCrsMat);
378  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
379 
383 
385  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xTpetraCrsMat);
386  return ret;
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 and LO=GO=int");
393  }
394 #endif
395  return Teuchos::null;
396  }
397 
399  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar> >& op) {
400 
401 #ifdef HAVE_XPETRA_TPETRA
402  if(isTpetra(op)) {
403  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
404  Teuchos::RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getTpetraOperator(op);
406  Teuchos::RCP<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
408  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
410 
414  return xTpetraCrsMat;
415  }
416 #endif
417 
418 #ifdef HAVE_XPETRA_EPETRA
419  if(isEpetra(op)) {
420  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double and LO=GO=int");
421  }
422 #endif
423  return Teuchos::null;
424  }
425 
428  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > thyraMap = Teuchos::null;
429 #ifdef HAVE_XPETRA_TPETRA
430  if(map->lib() == Xpetra::UseTpetra) {
432  if (tpetraMap == Teuchos::null)
433  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
434  RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetraMap->getTpetra_Map());
435  thyraMap = thyraTpetraMap;
436  }
437 #endif
438 
439 #ifdef HAVE_XPETRA_EPETRA
440  if(map->lib() == Xpetra::UseEpetra) {
441  TEUCHOS_TEST_FOR_EXCEPTION(thyraMap==Teuchos::null, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double and LO=GO=int");
442  }
443 #endif
444 
445  return thyraMap;
446  }
447 
448  // update Thyra multi vector with data from Xpetra multi vector
449  // In case of a Thyra::ProductMultiVector the Xpetra::MultiVector is splitted into its subparts using a provided MapExtractor
450  static void
451  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) {
452  using Teuchos::RCP;
453  using Teuchos::rcp_dynamic_cast;
454  using Teuchos::as;
455  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
456  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
457  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
458  //typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
459  //typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
460  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
462 
463  // copy data from tY_inout to Y_inout
464  RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
465  if(prodTarget != Teuchos::null) {
466  // SPECIAL CASE: product vector:
467  // update Thyra product multi vector with data from a merged Xpetra multi vector
468 
469  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
470  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.");
471 
472  for(int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
473  // access Xpetra data
474  RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
475 
476  // access Thyra data
477  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
478  RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
479  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
480  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
481  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
482  RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
483  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
484 
485  // loop over all vectors in multivector
486  for(size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
487  Teuchos::ArrayRCP< const Scalar > xpData = xpSubBlock->getData(j); // access const data from Xpetra object
488 
489  // loop over all local rows
490  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
491  (*thyData)(i,j) = xpData[i];
492  }
493  }
494  }
495  } else {
496  // STANDARD case:
497  // update Thyra::MultiVector with data from an Xpetra::MultiVector
498 
499  // access Thyra data
500  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
501  TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
502  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
503  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
504  RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
505  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
506 
507  // loop over all vectors in multivector
508  for(size_t j = 0; j < source->getNumVectors(); ++j) {
509  Teuchos::ArrayRCP< const Scalar > xpData = source->getData(j); // access const data from Xpetra object
510  // loop over all local rows
511  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
512  (*thyData)(i,j) = xpData[i];
513  }
514  }
515  }
516  }
517 
520  // create a Thyra operator from Xpetra::CrsMatrix
521  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
522 
523  bool bIsTpetra = false;
524 
525 #ifdef HAVE_XPETRA_TPETRA
527  if(tpetraMat!=Teuchos::null) {
528  bIsTpetra = true;
531  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
533 
534  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
536  Teuchos::RCP<const Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
538 
539  thyraOp = Thyra::createConstLinearOp(tpOperator);
541  }
542 #endif
543 
544 #ifdef HAVE_XPETRA_EPETRA
545  TEUCHOS_TEST_FOR_EXCEPTION(bIsTpetra == false, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double and LO=GO=int");
546 #endif
547  return thyraOp;
548  }
549 
552  // create a Thyra operator from Xpetra::CrsMatrix
553  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
554 
555  bool bIsTpetra = false;
556 
557 #ifdef HAVE_XPETRA_TPETRA
559  if(tpetraMat!=Teuchos::null) {
560  bIsTpetra = true;
563  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
565 
566  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
568  Teuchos::RCP<Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
570 
571  thyraOp = Thyra::createLinearOp(tpOperator);
573  }
574 #endif
575 
576 #ifdef HAVE_XPETRA_EPETRA
577  TEUCHOS_TEST_FOR_EXCEPTION(bIsTpetra == false, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double and LO=GO=int");
578 #endif
579  return thyraOp;
580  }
581 
584 
585  int nRows = mat->Rows();
586  int nCols = mat->Cols();
587 
589 
590  bool bTpetra = false;
591 #ifdef HAVE_XPETRA_TPETRA
593  if(tpetraMat!=Teuchos::null) bTpetra = true;
594 #endif
595 
596 #ifdef HAVE_XPETRA_EPETRA
597  TEUCHOS_TEST_FOR_EXCEPTION(bTpetra == false, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double and LO=GO=int");
598 #endif
599 
600  // create new Thyra blocked operator
602  Thyra::defaultBlockedLinearOp<Scalar>();
603 
604  blockMat->beginBlockFill(nRows,nCols);
605 
606  for (int r=0; r<nRows; ++r) {
607  for (int c=0; c<nCols; ++c) {
609  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(mat->getMatrix(r,c));
610  std::stringstream label; label << "A" << r << c;
611  blockMat->setBlock(r,c,thBlock);
612  }
613  }
614 
615  blockMat->endBlockFill();
616 
617  return blockMat;
618  }
619 
620 }; // end Utils class
621 
622 
623 // full specialization for Epetra support
624 // Note, that Thyra only has support for Epetra (not for Epetra64)
625 #ifdef HAVE_XPETRA_EPETRA
626 template <>
627 class ThyraUtils<double, int, int, EpetraNode> {
628 public:
629  typedef double Scalar;
630  typedef int LocalOrdinal;
631  typedef int GlobalOrdinal;
632  typedef EpetraNode Node;
633 
634 private:
635 #undef XPETRA_THYRAUTILS_SHORT
636 #include "Xpetra_UseShortNames.hpp"
637 
638 public:
639 
641  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) {
642 
644 
645  if(stridedBlockId == -1) {
646  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo.size() != 0);
647  }
648  else {
649  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo[stridedBlockId] != 0);
650  }
651 
653  return ret;
654  }
655 
657  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
658  using Teuchos::RCP;
659  using Teuchos::rcp_dynamic_cast;
660  using Teuchos::as;
661  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
662  typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
665  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
666 
667  RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase >(vectorSpace);
668  if(prodVectorSpace != Teuchos::null) {
669  // SPECIAL CASE: product Vector space
670  // merge all submaps to one large Xpetra Map
671  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error, "Found a product vector space with zero blocks.");
672  std::vector<RCP<Map> > maps(prodVectorSpace->numBlocks(), Teuchos::null);
673  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
674  RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
675  RCP<Map> map = ThyUtils::toXpetra(bv, comm); // recursive call
676  maps[b] = map;
677  }
678 
679  // get offsets for submap GIDs
680  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
681  for(int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
682  gidOffsets[i] = maps[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
683  }
684 
685  // loop over all sub maps and collect GIDs
686  std::vector<GlobalOrdinal> gids;
687  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b){
688  for(LocalOrdinal l = 0; l < as<LocalOrdinal>(maps[b]->getNodeNumElements()); ++l) {
689  GlobalOrdinal gid = maps[b]->getGlobalElement(l) + gidOffsets[b];
690  gids.push_back(gid);
691  }
692  }
693 
694  // create full map
696  Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
697  RCP<Map> fullMap = MapFactory::Build(maps[0]->lib(), INVALID, gidsView, maps[0]->getIndexBase(), comm);
698 
700  return fullMap;
701  } else {
702  // STANDARD CASE: no product map
703  // Epetra/Tpetra specific code to access the underlying map data
704 
705  // check whether we have a Tpetra based Thyra operator
706  bool bIsTpetra = false;
707 #ifdef HAVE_XPETRA_TPETRA
708 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
709  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
710  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vectorSpace);
711  bIsTpetra = Teuchos::is_null(tpetra_vsc) ? false : true;
712 #endif
713 #endif
714 
715  // check whether we have an Epetra based Thyra operator
716  bool bIsEpetra = !bIsTpetra; // note: this is a little bit fragile!
717 
718 #ifdef HAVE_XPETRA_TPETRA
719  if(bIsTpetra) {
720 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
721  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
722  typedef Thyra::VectorBase<Scalar> ThyVecBase;
723  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
724  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
725  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
726  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
728  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
730  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
732 
733  RCP<Map> rgXpetraMap = Xpetra::toXpetraNonConst(rgTpetraMap);
735  return rgXpetraMap;
736 #else
737  throw Xpetra::Exceptions::RuntimeError("Problem AAA. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
738 #endif
739  }
740 #endif
741 
742 #ifdef HAVE_XPETRA_EPETRA
743  if(bIsEpetra) {
744  // TODO check for product vector space!
745  RCP<const Epetra_Map> epMap = Teuchos::null;
746  RCP<const Epetra_Map>
747  epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map> >(vectorSpace,"epetra_map");
748  if(!Teuchos::is_null(epetra_map)) {
751  return rgXpetraMap;
752  } else {
753  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "No Epetra_Map data found in Thyra::VectorSpace.");
754  }
755  }
756 #endif
757  } // end standard case (no product map)
758  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "Cannot transform Thyra::VectorSpace to Xpetra::Map.");
759  return Teuchos::null;
760  }
761 
762  // const version
764  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
765  using Teuchos::RCP;
766  using Teuchos::rcp_dynamic_cast;
767  using Teuchos::as;
768  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
769  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
770  //typedef Thyra::VectorBase<Scalar> ThyVecBase;
771  //typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
772  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
774  //typedef Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node> MapFactory;
777  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
778 
779  // return value
780  RCP<MultiVector> xpMultVec = Teuchos::null;
781 
782  // check whether v is a product multi vector
783  Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase >(v);
784  if(thyProdVec != Teuchos::null) {
785  // SPECIAL CASE: product Vector
786  // merge all subvectors to one large Xpetra vector
787  RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
788 
789  // create new Epetra_MultiVector living on full map (stored in eMap)
790  xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
791 
792  // fill xpMultVec with Thyra data
793  std::vector<GlobalOrdinal> lidOffsets(thyProdVec->productSpace()->numBlocks()+1,0);
794  for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
795  Teuchos::RCP<const ThyVecSpaceBase> thySubMap = thyProdVec->productSpace()->getBlock(b);
796  Teuchos::RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(thySubMap);
797  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
798  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : thySubMap->dim() );
799  lidOffsets[b+1] = localSubDim + lidOffsets[b]; // calculate lid offset for next block
800  RCP<Thyra::ConstDetachedMultiVectorView<Scalar> > thyData =
801  Teuchos::rcp(new Thyra::ConstDetachedMultiVectorView<Scalar>(thyProdVec->getMultiVectorBlock(b),Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
802  for(size_t vv = 0; vv < xpMultVec->getNumVectors(); ++vv) {
803  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
804  xpMultVec->replaceLocalValue(i + lidOffsets[b] , vv, (*thyData)(i,vv));
805  }
806  }
807  }
809  return xpMultVec;
810  } else {
811  // STANDARD CASE: no product vector
812  // Epetra/Tpetra specific code to access the underlying map data
813  bool bIsTpetra = false;
814 #ifdef HAVE_XPETRA_TPETRA
815 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
816  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
817 
818  //typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
819  //typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
820  typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
821  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
822  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
824  typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
825 
826  RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
827  RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
828  if(thyraTpetraMultiVector != Teuchos::null) {
829  bIsTpetra = true;
830  const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
832  RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
833  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
834  xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
836  return xpMultVec;
837  }
838 #endif
839 #endif
840 
841 #ifdef HAVE_XPETRA_EPETRA
842  if(bIsTpetra == false) {
843  // no product vector
844  Teuchos::RCP<Map> map = ThyUtils::toXpetra(v->range(), comm);
846  RCP<Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = rcp_dynamic_cast<Xpetra::EpetraMapT<GlobalOrdinal,Node> >(map);
848  RCP<const Epetra_Map> eMap = Teuchos::rcpFromRef(xeMap->getEpetra_Map());
850  Teuchos::RCP<const Epetra_MultiVector> epMultVec = Thyra::get_Epetra_MultiVector(*eMap, v);
852  RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<Epetra_MultiVector>(epMultVec);
853  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epNonConstMultVec));
854  xpMultVec = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node>(epNonConstMultVec));
856  return xpMultVec;
857  }
858 #endif
859  } // end standard case
860  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "Cannot transform Thyra::MultiVector to Xpetra::MultiVector.");
861  return Teuchos::null;
862  }
863 
864  // non-const version
866  toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > v, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
868  Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar> >(v);
870  toXpetra(cv,comm);
871  return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(r);
872  }
873 
874  static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
875  // check whether we have a Tpetra based Thyra operator
876  bool bIsTpetra = false;
877 #ifdef HAVE_XPETRA_TPETRA
878 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
879  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
880 
881  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(op);
882  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
883 
884  // for debugging purposes: find out why dynamic cast failed
885  if(!bIsTpetra &&
886 #ifdef HAVE_XPETRA_EPETRA
887  Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
888 #endif
889  Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
890  // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
891  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
892  if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
893  std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
894  std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
895  std::cout << " properly set!" << std::endl;
896  std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
897  }
898  }
899 #endif
900 #endif
901 
902 #if 0
903  // Check whether it is a blocked operator.
904  // If yes, grab the (0,0) block and check the underlying linear algebra
905  // Note: we expect that the (0,0) block exists!
906  if(bIsTpetra == false) {
908  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
909  if(ThyBlockedOp != Teuchos::null) {
910  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
912  ThyBlockedOp->getBlock(0,0);
914  bIsTpetra = isTpetra(b00);
915  }
916  }
917 #endif
918 
919  return bIsTpetra;
920  }
921 
922  static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
923  // check whether we have an Epetra based Thyra operator
924  bool bIsEpetra = false;
925 
926 #ifdef HAVE_XPETRA_EPETRA
927  Teuchos::RCP<const Thyra::EpetraLinearOp> epetra_op = Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op,false);
928  bIsEpetra = Teuchos::is_null(epetra_op) ? false : true;
929 #endif
930 
931 #if 0
932  // Check whether it is a blocked operator.
933  // If yes, grab the (0,0) block and check the underlying linear algebra
934  // Note: we expect that the (0,0) block exists!
935  if(bIsEpetra == false) {
937  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op,false);
938  if(ThyBlockedOp != Teuchos::null) {
939  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
941  ThyBlockedOp->getBlock(0,0);
943  bIsEpetra = isEpetra(b00);
944  }
945  }
946 #endif
947 
948  return bIsEpetra;
949  }
950 
951  static bool isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
952  // Check whether it is a blocked operator.
954  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
955  if(ThyBlockedOp != Teuchos::null) {
956  return true;
957  }
958  return false;
959  }
960 
962  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >& op) {
963 
964 #ifdef HAVE_XPETRA_TPETRA
965  if(isTpetra(op)) {
966 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
967  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
968 
969  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
970  Teuchos::RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getConstTpetraOperator(op);
971  // we should also add support for the const versions!
972  //getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
974  Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
976  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
978  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraCrsMat);
979  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
980 
985  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xTpetraCrsMat);
987  return ret;
988 #else
989  throw Xpetra::Exceptions::RuntimeError("Problem BBB. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
990 #endif
991  }
992 #endif
993 
994 #ifdef HAVE_XPETRA_EPETRA
995  if(isEpetra(op)) {
996  Teuchos::RCP<const Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator( *op );
998  Teuchos::RCP<const Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_op);
1000  Teuchos::RCP<const Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(epetra_rowmat);
1001  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_crsmat));
1002  Teuchos::RCP<Epetra_CrsMatrix> epetra_ncnstcrsmat = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(epetra_crsmat);
1003  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_ncnstcrsmat));
1004 
1007  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1008 
1010  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xEpetraCrsMat);
1012  return ret;
1013  }
1014 #endif
1015  return Teuchos::null;
1016  }
1017 
1019  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar> >& op) {
1020 
1021 #ifdef HAVE_XPETRA_TPETRA
1022  if(isTpetra(op)) {
1023 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1024  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1025 
1026  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1027  Teuchos::RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getTpetraOperator(op);
1029  Teuchos::RCP<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
1031  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
1033 
1036  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
1037  return xTpetraCrsMat;
1038 #else
1039  throw Xpetra::Exceptions::RuntimeError("Problem CCC. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1040 #endif
1041  }
1042 #endif
1043 
1044 #ifdef HAVE_XPETRA_EPETRA
1045  if(isEpetra(op)) {
1046  Teuchos::RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator( *op );
1048  Teuchos::RCP<Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetra_op);
1049  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_rowmat));
1050  Teuchos::RCP<Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetra_rowmat);
1051  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_crsmat));
1052 
1055  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
1056  return xEpetraCrsMat;
1057  }
1058 #endif
1059  return Teuchos::null;
1060  }
1061 
1064  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > thyraMap = Teuchos::null;
1065 #ifdef HAVE_XPETRA_TPETRA
1066  if(map->lib() == Xpetra::UseTpetra) {
1067 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1068  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1070  if (tpetraMap == Teuchos::null)
1071  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
1072  RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > tpMap = tpetraMap->getTpetra_Map();
1073  RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
1074  thyraMap = thyraTpetraMap;
1075 #else
1076  throw Xpetra::Exceptions::RuntimeError("Problem DDD. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1077 #endif
1078  }
1079 #endif
1080 
1081 #ifdef HAVE_XPETRA_EPETRA
1082  if(map->lib() == Xpetra::UseEpetra) {
1084  if (epetraMap == Teuchos::null)
1085  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
1086  RCP<const Thyra::VectorSpaceBase<Scalar> > thyraEpetraMap = Thyra::create_VectorSpace(Teuchos::rcpFromRef(epetraMap->getEpetra_Map()));
1087  thyraMap = thyraEpetraMap;
1088  }
1089 #endif
1090 
1091  return thyraMap;
1092  }
1093 
1094  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) {
1095  using Teuchos::RCP;
1096  using Teuchos::rcp_dynamic_cast;
1097  using Teuchos::as;
1098  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1099  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
1100  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1101  //typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
1102  //typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
1103  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1105 
1106  // copy data from tY_inout to Y_inout
1107  RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
1108  if(prodTarget != Teuchos::null) {
1109  // SPECIAL CASE: product vector:
1110  // update Thyra product multi vector with data from a merged Xpetra multi vector
1111 
1112  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
1113  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.");
1114 
1115  for(int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1116  // access Xpetra data
1117  RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
1118 
1119  // access Thyra data
1120  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
1121  RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
1122  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
1123  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
1124  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
1125  RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1126  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
1127 
1128  // loop over all vectors in multivector
1129  for(size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
1130  Teuchos::ArrayRCP< const Scalar > xpData = xpSubBlock->getData(j); // access const data from Xpetra object
1131 
1132  // loop over all local rows
1133  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1134  (*thyData)(i,j) = xpData[i];
1135  }
1136  }
1137  }
1138  } else {
1139  // STANDARD case:
1140  // update Thyra::MultiVector with data from an Xpetra::MultiVector
1141 
1142  // access Thyra data
1143  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
1144  TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
1145  const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
1146  const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
1147  RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1148  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target,Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
1149 
1150  // loop over all vectors in multivector
1151  for(size_t j = 0; j < source->getNumVectors(); ++j) {
1152  Teuchos::ArrayRCP< const Scalar > xpData = source->getData(j); // access const data from Xpetra object
1153  // loop over all local rows
1154  for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1155  (*thyData)(i,j) = xpData[i];
1156  }
1157  }
1158  }
1159  }
1160 
1163  // create a Thyra operator from Xpetra::CrsMatrix
1164  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
1165 
1166 #ifdef HAVE_XPETRA_TPETRA
1168  if(tpetraMat!=Teuchos::null) {
1169 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1170  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1171 
1174  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
1176 
1177  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
1179  Teuchos::RCP<const Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
1181 
1182  thyraOp = Thyra::createConstLinearOp(tpOperator);
1184 #else
1185  throw Xpetra::Exceptions::RuntimeError("Problem EEE. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1186 #endif
1187  }
1188 #endif
1189 
1190 #ifdef HAVE_XPETRA_EPETRA
1192  if(epetraMat!=Teuchos::null) {
1195  Teuchos::RCP<const Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrix();
1197 
1198  Teuchos::RCP<const Thyra::EpetraLinearOp> thyraEpOp = Thyra::epetraLinearOp(epCrsMat,"op");
1200  thyraOp = thyraEpOp;
1201  }
1202 #endif
1203  return thyraOp;
1204  }
1205 
1208  // create a Thyra operator from Xpetra::CrsMatrix
1209  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
1210 
1211 #ifdef HAVE_XPETRA_TPETRA
1213  if(tpetraMat!=Teuchos::null) {
1214 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1215  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1216 
1219  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
1221 
1222  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
1224  Teuchos::RCP<Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
1226 
1227  thyraOp = Thyra::createLinearOp(tpOperator);
1229 #else
1230  throw Xpetra::Exceptions::RuntimeError("Problem FFF. Add TPETRA_INST_INT_INT:BOOL=ON in your configuration.");
1231 #endif
1232  }
1233 #endif
1234 
1235 #ifdef HAVE_XPETRA_EPETRA
1237  if(epetraMat!=Teuchos::null) {
1240  Teuchos::RCP<Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrixNonConst();
1242 
1243  Teuchos::RCP<Thyra::EpetraLinearOp> thyraEpOp = Thyra::nonconstEpetraLinearOp(epCrsMat,"op");
1245  thyraOp = thyraEpOp;
1246  }
1247 #endif
1248  return thyraOp;
1249  }
1250 
1253 #if 0
1254  {
1255 
1256  int nRows = mat->Rows();
1257  int nCols = mat->Cols();
1258 
1260 
1261  bool bTpetra = false;
1262  bool bEpetra = false;
1263 #ifdef HAVE_XPETRA_TPETRA
1264 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1265  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1267  if(tpetraMat!=Teuchos::null) bTpetra = true;
1268 #else
1269  bTpetra = false;
1270 #endif
1271 #endif
1272 
1273 #ifdef HAVE_XPETRA_EPETRA
1275  if(epetraMat!=Teuchos::null) bEpetra = true;
1276 #endif
1277 
1278  TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra
1279 
1280  // create new Thyra blocked operator
1282  Thyra::defaultBlockedLinearOp<Scalar>();
1283 
1284  blockMat->beginBlockFill(nRows,nCols);
1285 
1286  for (int r=0; r<nRows; ++r) {
1287  for (int c=0; c<nCols; ++c) {
1289  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(mat->getMatrix(r,c));
1290  std::stringstream label; label << "A" << r << c;
1291  blockMat->setBlock(r,c,thBlock);
1292  }
1293  }
1294 
1295  blockMat->endBlockFill();
1296 
1297  return blockMat;
1298  }
1299 #endif
1300 
1301 }; // specialization on SC=double, LO=GO=int and NO=EpetraNode
1302 #endif // HAVE_XPETRA_EPETRA
1303 } // end namespace Xpetra
1304 
1305 #define XPETRA_THYRAUTILS_SHORT
1306 #endif // HAVE_XPETRA_THYRA
1307 
1308 #endif // XPETRA_THYRAUTILS_HPP
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
GlobalOrdinal GO
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Xpetra namespace
Exception throws to report errors in the internal logical of the program.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
Definition: Xpetra_Map.hpp:68
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static RCP< StridedMap > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
const RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > toXpetraNonConst(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &)
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
TypeTo as(const TypeFrom &t)
Create an Xpetra::Map instance.
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.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)