Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Xpetra_ThyraUtils_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Xpetra: A linear algebra interface package
4 //
5 // Copyright 2012 NTESS and the Xpetra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef HAVE_XPETRA_THYRAUTILS_DEF_HPP
11 #define HAVE_XPETRA_THYRAUTILS_DEF_HPP
12 
13 #ifdef HAVE_XPETRA_THYRA
14 
15 #include "Xpetra_BlockedCrsMatrix.hpp"
16 
17 #include "Xpetra_ThyraUtils_decl.hpp"
18 
19 namespace Xpetra {
20 
21 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
22 Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>
23 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
24  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm, std::vector<size_t>& stridingInfo, LocalOrdinal stridedBlockId, GlobalOrdinal offset) {
25  Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map = toXpetra(vectorSpace, comm);
26 
27  if (stridedBlockId == -1) {
28  TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo.size() != 0);
29  } else {
30  TEUCHOS_TEST_FOR_EXCEPT(map->getLocalNumElements() % stridingInfo[stridedBlockId] != 0);
31  }
32 
33  Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>> ret = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(map, stridingInfo, stridedBlockId, offset);
34  return ret;
35 }
36 
37 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
38 Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
39 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
40  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
41  using Teuchos::as;
42  using Teuchos::RCP;
43  using Teuchos::rcp_dynamic_cast;
44  using ThyVecSpaceBase = Thyra::VectorSpaceBase<Scalar>;
45  using ThyProdVecSpaceBase = Thyra::ProductVectorSpaceBase<Scalar>;
46  using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
47 
48  RCP<Map> resultMap = Teuchos::null;
49  RCP<const ThyProdVecSpaceBase> prodVectorSpace = rcp_dynamic_cast<const ThyProdVecSpaceBase>(vectorSpace);
50  if (prodVectorSpace != Teuchos::null) {
51  // SPECIAL CASE: product Vector space
52  // collect all submaps to store them in a hierarchical BlockedMap object
53  TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks() == 0, std::logic_error, "Found a product vector space with zero blocks.");
54  std::vector<RCP<const Map>> mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
55  std::vector<RCP<const Map>> mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
56  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
57  RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
58  // can be of type Map or BlockedMap (containing Thyra GIDs)
59  mapsThyra[b] = ThyUtils::toXpetra(bv, comm); // recursive call
60  }
61 
62  // get offsets for submap GIDs
63  // we need that for the full map (Xpetra GIDs)
64  std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(), 0);
65  for (int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
66  gidOffsets[i] = mapsThyra[i - 1]->getMaxAllGlobalIndex() + gidOffsets[i - 1] + 1;
67  }
68 
69  for (int b = 0; b < prodVectorSpace->numBlocks(); ++b) {
70  RCP<const ThyVecSpaceBase> bv = prodVectorSpace->getBlock(b);
71  // map can be of type Map or BlockedMap (containing Xpetra style GIDs)
72  mapsXpetra[b] = MapUtils::transformThyra2XpetraGIDs(*mapsThyra[b], gidOffsets[b]);
73  }
74 
75  resultMap = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>(mapsXpetra, mapsThyra));
76  } else {
77 #ifdef HAVE_XPETRA_TPETRA
78  // STANDARD CASE: no product map
79  // check whether we have a Tpetra based Thyra operator
80  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(vectorSpace);
81  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.");
82 
83  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> TpMap;
84  typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpVector;
85  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
86  typedef Thyra::VectorBase<Scalar> ThyVecBase;
87  RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
88  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
89  RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
90  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
91  RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
92  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
93  resultMap = Xpetra::toXpetraNonConst(rgTpetraMap);
94 #else
95  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.");
96 #endif
97  } // end standard case (no product map)
98  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(resultMap));
99  return resultMap;
100 }
101 
102 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
103 Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
104 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
105  toXpetra(Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
106  using Teuchos::as;
107  using Teuchos::RCP;
108  using Teuchos::rcp_dynamic_cast;
109 
110  using ThyMultVecBase = Thyra::MultiVectorBase<Scalar>;
111  using ThyProdMultVecBase = Thyra::ProductMultiVectorBase<Scalar>;
112  using ThyUtils = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
113 
114  // return value
115  RCP<MultiVector> xpMultVec = Teuchos::null;
116 
117  // check whether v is a product multi vector
118  Teuchos::RCP<const ThyProdMultVecBase> thyProdVec = rcp_dynamic_cast<const ThyProdMultVecBase>(v);
119  if (thyProdVec != Teuchos::null) {
120  // SPECIAL CASE: create a nested BlockedMultiVector
121  // generate nested BlockedMap (containing Thyra and Xpetra GIDs)
122  RCP<Map> fullMap = ThyUtils::toXpetra(v->range(), comm);
123  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(Teuchos::rcp_dynamic_cast<BlockedMap>(fullMap)));
124 
125  // create new Xpetra::BlockedMultiVector
126  xpMultVec = MultiVectorFactory::Build(fullMap, as<size_t>(thyProdVec->domain()->dim()));
127 
128  RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec, true);
129 
130  // loop over all blocks, transform Thyra MultiVectors to Xpetra MultiVectors recursively
131  for (int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b) {
132  RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
133  // xpBlockMV can be of type MultiVector or BlockedMultiVector
134  RCP<const MultiVector> xpBlockMV = ThyUtils::toXpetra(thyBlockMV, comm); // recursive call
135  xpBlockedMultVec->setMultiVector(b, xpBlockMV, true /* Thyra mode */);
136  }
137  } else {
138  // STANDARD CASE: no product vector
139 #ifdef HAVE_XPETRA_TPETRA
140  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> ConverterT;
141  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpMultVec;
142  typedef Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpTpMultVec;
143  typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
144  typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
145 
146  RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<const ThySpmdMultVecBase>(v);
147  RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<const ThyTpMultVec>(thyraSpmdMultiVector);
148  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.");
149  const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
150  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpMultVec));
151  RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
152  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpNonConstMultVec));
153  xpMultVec = rcp(new XpTpMultVec(tpNonConstMultVec));
154 #else
155  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.");
156 #endif
157  } // end standard case
158  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpMultVec));
159  return xpMultVec;
160 }
161 
162 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
163 Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
164 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
165  toXpetra(Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> v, const Teuchos::RCP<const Teuchos::Comm<int>>& comm) {
166  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> cv =
167  Teuchos::rcp_const_cast<const Thyra::MultiVectorBase<Scalar>>(v);
168  Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> r =
169  toXpetra(cv, comm);
170  return Teuchos::rcp_const_cast<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(r);
171 }
172 
173 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
174 bool Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
175  isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
176  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(op));
177 
178  // check whether we have a Tpetra based Thyra operator
179  bool bIsTpetra = false;
180 #ifdef HAVE_XPETRA_TPETRA
181  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(op);
182  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
183 
184  // for debugging purposes: find out why dynamic cast failed
185  if (!bIsTpetra &&
186 #ifdef HAVE_XPETRA_EPETRA
187  Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
188 #endif
189  Teuchos::rcp_dynamic_cast<const Thyra::DefaultBlockedLinearOp<Scalar>>(op) == Teuchos::null) {
190  // If op is not blocked and not an Epetra object, it should be in fact an Tpetra object
191  typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
192  if (Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
193  std::cout << "ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
194  std::cout << " If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
195  std::cout << " properly set!" << std::endl;
196  std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op, true) << std::endl;
197  }
198  }
199 #endif
200 
201 #if 0
202  // Check whether it is a blocked operator.
203  // If yes, grab the (0,0) block and check the underlying linear algebra
204  // Note: we expect that the (0,0) block exists!
205  if(bIsTpetra == false) {
206  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
207  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(op);
208  if(ThyBlockedOp != Teuchos::null) {
209  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
210  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > b00 =
211  ThyBlockedOp->getBlock(0,0);
212  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
213  bIsTpetra = isTpetra(b00);
214  }
215  }
216 #endif
217 
218  return bIsTpetra;
219 }
220 
221 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
222 bool Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
223  isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
224  return false;
225 }
226 
227 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
228 bool Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
229  isBlockedOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
230  // Check whether it is a blocked operator.
231  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar>> ThyBlockedOp =
232  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar>>(op);
233  if (ThyBlockedOp != Teuchos::null) {
234  return true;
235  }
236  return false;
237 }
238 
239 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
240 Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
241 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
242  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
243 #ifdef HAVE_XPETRA_TPETRA
244  if (isTpetra(op)) {
245  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
246  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
247  // we should also add support for the const versions!
248  // getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
249  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
250  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
251  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
252  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat);
253  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
254 
255  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
256  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraNcnstCrsMat));
257  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
258 
259  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
260  Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
261  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
262  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
263  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
264  Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
265  return xpMat;
266  }
267 #endif
268 
269 #ifdef HAVE_XPETRA_EPETRA
270  if (isEpetra(op)) {
271  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
272  }
273 #endif
274  return Teuchos::null;
275 }
276 
277 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
278 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
279 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
280  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
281  using Teuchos::RCP;
282  using Teuchos::rcp;
283  using Teuchos::rcp_const_cast;
284  using Teuchos::rcp_dynamic_cast;
285 
286 #ifdef HAVE_XPETRA_TPETRA
287  if (isTpetra(op)) {
288  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
289  typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraOperator_t;
290  RCP<const TpetraOperator_t> TpetraOp = TOE::getConstTpetraOperator(op);
291  TEUCHOS_TEST_FOR_EXCEPT(is_null(TpetraOp));
292  RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraRowMat =
293  rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp, true);
294  RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraCrsMat =
295  rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraRowMat, true);
296 
297  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraCrsMat =
298  rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(
299  rcp_const_cast<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraCrsMat)));
300  TEUCHOS_TEST_FOR_EXCEPT(is_null(xTpetraCrsMat));
301  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsMat =
302  rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraCrsMat, true);
303  RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpCrsWrap =
304  rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrsMat));
305  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpMat =
306  rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpCrsWrap, true);
307  return xpMat;
308  }
309 #endif
310 
311 #ifdef HAVE_XPETRA_EPETRA
312  if (isEpetra(op)) {
313  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
314  }
315 #endif
316  return Teuchos::null;
317 }
318 
319 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
320 Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
321 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
322  toXpetraOperator(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>& op) {
323 #ifdef HAVE_XPETRA_TPETRA
324  if (isTpetra(op)) {
325  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
326  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getConstTpetraOperator(op);
327 
328  auto nonConstTpetraOp = Teuchos::rcp_const_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(TpetraOp);
329 
330  Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
331  Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(nonConstTpetraOp));
332  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
333 
334  Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
335  Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
336  return xpOp;
337  }
338 #endif
339 
340 #ifdef HAVE_XPETRA_EPETRA
341  if (isEpetra(op)) {
342  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
343  }
344 #endif
345  return Teuchos::null;
346 }
347 
348 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
349 Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
350 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
351  toXpetraOperator(const Teuchos::RCP<Thyra::LinearOpBase<Scalar>>& op) {
352 #ifdef HAVE_XPETRA_TPETRA
353  if (isTpetra(op)) {
354  typedef Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node> TOE;
355  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> TpetraOp = TOE::getTpetraOperator(op);
356 
357  Teuchos::RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpetraOp =
358  Teuchos::rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(TpetraOp));
359  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraOp));
360 
361  Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpOp =
362  Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xTpetraOp, true);
363  return xpOp;
364  }
365 #endif
366 
367 #ifdef HAVE_XPETRA_EPETRA
368  if (isEpetra(op)) {
369  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
370  }
371 #endif
372  return Teuchos::null;
373 }
374 
375 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
376 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
377 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
378  toXpetra(const Teuchos::RCP<Thyra::DiagonalLinearOpBase<Scalar>>& op) {
379  using Teuchos::rcp_const_cast;
380  using Teuchos::rcp_dynamic_cast;
381 
382  RCP<const Thyra::VectorBase<Scalar>> diag = op->getDiag();
383 
384  RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpDiag;
385 #ifdef HAVE_XPETRA_TPETRA
386  using thyTpV = Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
387  using tV = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
388  if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
389  RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
390  if (!tDiag.is_null())
391  xpDiag = Xpetra::toXpetra(tDiag);
392  }
393 #endif
394  TEUCHOS_ASSERT(!xpDiag.is_null());
395  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
396  return M;
397 }
398 
399 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
400 Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
401 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
402  toXpetra(const Teuchos::RCP<const Thyra::DiagonalLinearOpBase<Scalar>>& op) {
403  return toXpetra(Teuchos::rcp_const_cast<Thyra::DiagonalLinearOpBase<Scalar>>(op));
404 }
405 
406 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
407 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
408 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
409  toThyra(Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> map) {
410  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> thyraMap = Teuchos::null;
411 
412  // check whether map is of type BlockedMap
413  RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
414  if (bmap.is_null() == false) {
415  Teuchos::Array<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>> vecSpaces(bmap->getNumMaps());
416  for (size_t i = 0; i < bmap->getNumMaps(); i++) {
417  // we need Thyra GIDs for all the submaps
418  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>> vs =
419  Xpetra::ThyraUtils<Scalar, LO, GO, Node>::toThyra(bmap->getMap(i, true));
420  vecSpaces[i] = vs;
421  }
422 
423  thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
424  return thyraMap;
425  }
426 
427  // standard case
428 #ifdef HAVE_XPETRA_TPETRA
429  if (map->lib() == Xpetra::UseTpetra) {
430  Teuchos::RCP<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>> tpetraMap = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal, GlobalOrdinal, Node>>(map);
431  if (tpetraMap == Teuchos::null)
432  throw Exceptions::BadCast("Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
433  RCP<Thyra::TpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>> thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetraMap->getTpetra_Map());
434  thyraMap = thyraTpetraMap;
435  }
436 #endif
437 
438 #ifdef HAVE_XPETRA_EPETRA
439  if (map->lib() == Xpetra::UseEpetra) {
440  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
441  }
442 #endif
443 
444  return thyraMap;
445 }
446 
447 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
448 Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>>
449 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
450  toThyraMultiVector(Teuchos::RCP<const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
451  // create Thyra MultiVector
452 #ifdef HAVE_XPETRA_TPETRA
453  if (vec->getMap()->lib() == Xpetra::UseTpetra) {
454  auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
455  RCP<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpMV = Teuchos::rcp_dynamic_cast<const TpetraMultiVector>(vec)->getTpetra_MultiVector();
456  auto thyDomMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(vec->getNumVectors(), vec->getMap()->getComm()));
457  auto thyMV = rcp(new Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
458  thyMV->initialize(thyTpMap, thyDomMap, tpMV);
459  return thyMV;
460  }
461 #endif
462 
463 #ifdef HAVE_XPETRA_EPETRA
464  if (vec->getMap()->lib() == Xpetra::UseEpetra) {
465  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
466  }
467 #endif
468 
469  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "MultiVector cannot be converted to Thyra.");
470 }
471 
472 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
473 Teuchos::RCP<const Thyra::VectorBase<Scalar>>
474 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
475  toThyraVector(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> vec) {
476  // create Thyra Vector
477 #ifdef HAVE_XPETRA_TPETRA
478  if (vec->getMap()->lib() == Xpetra::UseTpetra) {
479  auto thyTpMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(Teuchos::rcp_dynamic_cast<const TpetraMap>(vec->getMap())->getTpetra_Map());
480  RCP<Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpVec = Teuchos::rcp_dynamic_cast<const TpetraVector>(vec)->getTpetra_Vector();
481  auto thyVec = rcp(new Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>());
482  thyVec->initialize(thyTpMap, tpVec);
483  return thyVec;
484  }
485 #endif
486 
487 #ifdef HAVE_XPETRA_EPETRA
488  if (vec->getMap()->lib() == Xpetra::UseEpetra) {
489  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
490  }
491 #endif
492 
493  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Vector cannot be converted to Thyra.");
494 }
495 
496 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
497 void Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
498  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) {
499  using Teuchos::as;
500  using Teuchos::RCP;
501  using Teuchos::rcp_dynamic_cast;
502  typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
503  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
504  typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
505  // typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
506  // typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
507  typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
508 
509  // copy data from tY_inout to Y_inout
510  RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
511  if (prodTarget != Teuchos::null) {
512  RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<const BlockedMultiVector>(source);
513  if (bSourceVec.is_null() == true) {
514  // SPECIAL CASE: target vector is product vector:
515  // update Thyra product multi vector with data from a merged Xpetra multi vector
516 
517  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error, "Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
518  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.");
519 
520  for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
521  // access Xpetra data
522  RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb, false); // use Xpetra ordering (doesn't really matter)
523 
524  // access Thyra data
525  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
526  RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
527  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(vs);
528  const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
529  const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim());
530  RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
531  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*thySubBlock, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
532 
533  // loop over all vectors in multivector
534  for (size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
535  Teuchos::ArrayRCP<const Scalar> xpData = xpSubBlock->getData(j); // access const data from Xpetra object
536 
537  // loop over all local rows
538  for (LocalOrdinal i = 0; i < localSubDim; ++i) {
539  (*thyData)(i, j) = xpData[i];
540  }
541  }
542  }
543  } else {
544  // source vector is a blocked multivector
545  // TODO test me
546  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.");
547 
548  for (int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
549  // access Thyra data
550  RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb, true); // use Thyra ordering
551 
552  Teuchos::RCP<const ThyMultVecBase> thyXpSubBlock = toThyraMultiVector(xpSubBlock);
553 
554  // access Thyra data
555  Teuchos::RCP<ThyMultVecBase> thySubBlock = prodTarget->getNonconstMultiVectorBlock(bbb);
556  Thyra::assign(thySubBlock.ptr(), *thyXpSubBlock);
557  }
558  }
559  } else {
560  // STANDARD case:
561  // update Thyra::MultiVector with data from an Xpetra::MultiVector
562 
563  // access Thyra data
564  RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<const ThySpmdVecSpaceBase>(target->range());
565  TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error, "Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
566  const LocalOrdinal localOffset = (mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0);
567  const LocalOrdinal localSubDim = (mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim());
568  RCP<Thyra::DetachedMultiVectorView<Scalar>> thyData =
569  Teuchos::rcp(new Thyra::DetachedMultiVectorView<Scalar>(*target, Teuchos::Range1D(localOffset, localOffset + localSubDim - 1)));
570 
571  // loop over all vectors in multivector
572  for (size_t j = 0; j < source->getNumVectors(); ++j) {
573  Teuchos::ArrayRCP<const Scalar> xpData = source->getData(j); // access const data from Xpetra object
574  // loop over all local rows
575  for (LocalOrdinal i = 0; i < localSubDim; ++i) {
576  (*thyData)(i, j) = xpData[i];
577  }
578  }
579  }
580 }
581 
582 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
583 Teuchos::RCP<const Thyra::LinearOpBase<Scalar>>
584 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
585  toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
586  // create a Thyra operator from Xpetra::CrsMatrix
587  Teuchos::RCP<const Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
588 
589  // bool bIsTpetra = false;
590 
591 #ifdef HAVE_XPETRA_TPETRA
592  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
593  if (tpetraMat != Teuchos::null) {
594  // bIsTpetra = true;
595  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
596  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
597  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
598 
599  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
600  Teuchos::RCP<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
601 
602  thyraOp = Thyra::createConstLinearOp(tpOperator);
603  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
604  } else {
605 #ifdef HAVE_XPETRA_EPETRA
606  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");
607 #else
608  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to Tpetra::CrsMatrix failed. Assume matrix should be Epetra then. No Epetra available");
609 #endif
610  }
611  return thyraOp;
612 #else
613  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
614  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
615 #endif
616 }
617 
618 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
619 Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
620 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
621  toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
622  // create a Thyra operator from Xpetra::CrsMatrix
623  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thyraOp = Teuchos::null;
624 
625  // bool bIsTpetra = false;
626 
627 #ifdef HAVE_XPETRA_TPETRA
628  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat);
629  if (tpetraMat != Teuchos::null) {
630  // bIsTpetra = true;
631  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(mat, true);
632  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
633  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
634 
635  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsMat, true);
636  Teuchos::RCP<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpRowMat, true);
637 
638  thyraOp = Thyra::createLinearOp(tpOperator);
639  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
640  } else {
641  // cast to TpetraCrsMatrix failed
642 #ifdef HAVE_XPETRA_EPETRA
643  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");
644 #else
645  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Guess, matrix should be Epetra then, but no Epetra available.");
646 #endif
647  }
648  return thyraOp;
649 #else
650  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
651  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
652 #endif
653 }
654 
655 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
656 Teuchos::RCP<Thyra::LinearOpBase<Scalar>>
657 Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
658  toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& mat) {
659  int nRows = mat->Rows();
660  int nCols = mat->Cols();
661 
662  Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Ablock = mat->getInnermostCrsMatrix();
663  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock);
664  TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
665 
666 #ifdef HAVE_XPETRA_TPETRA
667  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
668  if (tpetraMat != Teuchos::null) {
669  // create new Thyra blocked operator
670  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar>> blockMat =
671  Thyra::defaultBlockedLinearOp<Scalar>();
672 
673  blockMat->beginBlockFill(nRows, nCols);
674 
675  for (int r = 0; r < nRows; ++r) {
676  for (int c = 0; c < nCols; ++c) {
677  Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r, c);
678 
679  if (xpmat == Teuchos::null) continue; // shortcut for empty blocks
680 
681  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thBlock = Teuchos::null;
682 
683  // check whether the subblock is again a blocked operator
684  Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xpblock =
685  Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpmat);
686  if (xpblock != Teuchos::null) {
687  if (xpblock->Rows() == 1 && xpblock->Cols() == 1) {
688  // If it is a single block operator, unwrap it
689  Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
690  TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
691  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
692  } else {
693  // recursive call for general blocked operators
694  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpblock);
695  }
696  } else {
697  // check whether it is a CRSMatrix object
698  Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
699  TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
700  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
701  }
702 
703  blockMat->setBlock(r, c, thBlock);
704  }
705  }
706 
707  blockMat->endBlockFill();
708 
709  return blockMat;
710  } else {
711  // tpetraMat == Teuchos::null
712 #ifdef HAVE_XPETRA_EPETRA
713  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");
714 #else
715  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Cast to TpetraCrsMatrix failed. Guess, matrix should be Epetra then, but no Epetra available.");
716 #endif
717  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
718  }
719 #endif // endif HAVE_XPETRA_TPETRA
720 
721  // 4-Aug-2017 JJH Added 2nd condition to avoid "warning: dynamic initialization in unreachable code"
722  // If HAVE_XPETRA_TPETRA is defined, then this method will always return or throw in the if-then-else above.
723 #if defined(HAVE_XPETRA_EPETRA) && !defined(HAVE_XPETRA_TPETRA)
724  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Epetra needs SC=double, LO=int, and GO=int or GO=long long");
725  TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
726 #endif // endif HAVE_XPETRA_EPETRA
727 }
728 
729 #ifdef HAVE_XPETRA_EPETRA
730 
731 #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
732 // implementation of "toThyra" for full specialization on SC=double, LO=GO=int and NO=EpetraNode
733 // We need the specialization in the cpp file due to a circle dependency in the .hpp files for BlockedCrsMatrix
734 Teuchos::RCP<Thyra::LinearOpBase<double>>
735 ThyraUtils<double, int, int, EpetraNode>::toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, int, EpetraNode>>& mat) {
736  int nRows = mat->Rows();
737  int nCols = mat->Cols();
738 
739  Teuchos::RCP<Xpetra::Matrix<double, int, int, EpetraNode>> Ablock = mat->getInnermostCrsMatrix();
740  Teuchos::RCP<Xpetra::CrsMatrixWrap<double, int, int, EpetraNode>> Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<double, int, int, EpetraNode>>(Ablock);
741  TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
742 
743  bool bTpetra = false;
744  bool bEpetra = false;
745 #ifdef HAVE_XPETRA_TPETRA
746  // Note: Epetra is enabled
747 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
748  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
749  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
750  if (tpetraMat != Teuchos::null) bTpetra = true;
751 #else
752  bTpetra = false;
753 #endif
754 #endif
755 
756 #ifdef HAVE_XPETRA_EPETRA
757  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
758  if (epetraMat != Teuchos::null) bEpetra = true;
759 #endif
760 
761  TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra
762 
763  // create new Thyra blocked operator
764  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar>> blockMat =
765  Thyra::defaultBlockedLinearOp<Scalar>();
766 
767  blockMat->beginBlockFill(nRows, nCols);
768 
769  for (int r = 0; r < nRows; ++r) {
770  for (int c = 0; c < nCols; ++c) {
771  Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r, c);
772 
773  if (xpmat == Teuchos::null) continue; // shortcut for empty blocks
774 
775  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thBlock = Teuchos::null;
776 
777  // check whether the subblock is again a blocked operator
778  Teuchos::RCP<BlockedCrsMatrix> xpblock = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(xpmat);
779  if (xpblock != Teuchos::null) {
780  if (xpblock->Rows() == 1 && xpblock->Cols() == 1) {
781  // If it is a single block operator, unwrap it
782  Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
783  TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
784  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
785  } else {
786  // recursive call for general blocked operators
787  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpblock);
788  }
789  } else {
790  // check whether it is a CRSMatrix object
791  Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
792  TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
793  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
794  }
795 
796  blockMat->setBlock(r, c, thBlock);
797  }
798  }
799 
800  blockMat->endBlockFill();
801 
802  return blockMat;
803 }
804 #endif // #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
805 
806 #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
807 // implementation of "toThyra" for full specialization on SC=double, LO=int, GO=long long and NO=EpetraNode
808 // We need the specialization in the cpp file due to a circle dependency in the .hpp files for BlockedCrsMatrix
809 Teuchos::RCP<Thyra::LinearOpBase<double>>
810 ThyraUtils<double, int, long long, EpetraNode>::toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, long long, EpetraNode>>& mat) {
811  int nRows = mat->Rows();
812  int nCols = mat->Cols();
813 
814  Teuchos::RCP<Xpetra::Matrix<double, int, long long, EpetraNode>> Ablock = mat->getInnermostCrsMatrix();
815  Teuchos::RCP<Xpetra::CrsMatrixWrap<double, int, long long, EpetraNode>> Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<double, int, long long, EpetraNode>>(Ablock);
816  TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);
817 
818  bool bTpetra = false;
819  bool bEpetra = false;
820 #ifdef HAVE_XPETRA_TPETRA
821  // Note: Epetra is enabled
822 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
823  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
824  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
825  if (tpetraMat != Teuchos::null) bTpetra = true;
826 #else
827  bTpetra = false;
828 #endif
829 #endif
830 
831 #ifdef HAVE_XPETRA_EPETRA
832  Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(Ablock_wrap->getCrsMatrix());
833  if (epetraMat != Teuchos::null) bEpetra = true;
834 #endif
835 
836  TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra
837 
838  // create new Thyra blocked operator
839  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar>> blockMat =
840  Thyra::defaultBlockedLinearOp<Scalar>();
841 
842  blockMat->beginBlockFill(nRows, nCols);
843 
844  for (int r = 0; r < nRows; ++r) {
845  for (int c = 0; c < nCols; ++c) {
846  Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r, c);
847 
848  if (xpmat == Teuchos::null) continue; // shortcut for empty blocks
849 
850  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> thBlock = Teuchos::null;
851 
852  // check whether the subblock is again a blocked operator
853  Teuchos::RCP<BlockedCrsMatrix> xpblock = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(xpmat);
854  if (xpblock != Teuchos::null) {
855  if (xpblock->Rows() == 1 && xpblock->Cols() == 1) {
856  // If it is a single block operator, unwrap it
857  Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpblock->getCrsMatrix());
858  TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
859  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
860  } else {
861  // recursive call for general blocked operators
862  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpblock);
863  }
864  } else {
865  // check whether it is a CRSMatrix object
866  Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
867  TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);
868  thBlock = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpwrap->getCrsMatrix());
869  }
870 
871  blockMat->setBlock(r, c, thBlock);
872  }
873  }
874 
875  blockMat->endBlockFill();
876 
877  return blockMat;
878 }
879 #endif // #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
880 
881 #endif
882 
883 } // namespace Xpetra
884 
885 #endif
886 
887 #endif
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:24
Xpetra::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: Xpetra_Exceptions.hpp:65
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:191
Xpetra_ThyraUtils_decl.hpp
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
Definition: Xpetra_MatrixFactory_def.hpp:22
Xpetra::TpetraCrsMatrix::getTpetra_CrsMatrix
RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrix() const
Get the underlying Tpetra matrix.
Definition: Xpetra_TpetraCrsMatrix_def.hpp:611
Xpetra::TpetraCrsMatrix::getTpetra_CrsMatrixNonConst
RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetra_CrsMatrixNonConst() const
Get the underlying Tpetra matrix.
Definition: Xpetra_TpetraCrsMatrix_def.hpp:615
Xpetra::CrsMatrixWrap
Concrete implementation of Xpetra::Matrix.
Definition: Xpetra_CrsMatrixWrap_decl.hpp:46
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:24
Xpetra::Matrix
Xpetra-specific matrix class.
Definition: Xpetra_Matrix_decl.hpp:60
Xpetra::toXpetraNonConst
const RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > toXpetraNonConst(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &map)
Definition: Xpetra_TpetraMap_decl.hpp:228
Xpetra::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Definition: Xpetra_EpetraCrsGraph.cpp:29