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 #ifdef HAVE_XPETRA_THYRA
51 
52 #include "Xpetra_ConfigDefs.hpp"
53 
54 #ifdef HAVE_XPETRA_TPETRA
55 #include "Tpetra_ConfigDefs.hpp"
56 #endif
57 
58 #ifdef HAVE_XPETRA_EPETRA
59 #include "Epetra_CombineMode.h"
60 #endif
61 
62 #include "Xpetra_Map.hpp"
63 #include "Xpetra_StridedMap.hpp"
64 #include "Xpetra_StridedMapFactory.hpp"
65 #include "Xpetra_MapExtractor.hpp"
66 
67 #include <Thyra_VectorSpaceBase.hpp>
68 #include <Thyra_ProductVectorSpaceBase.hpp>
69 #include <Thyra_VectorSpaceBase.hpp>
70 #include <Thyra_DefaultBlockedLinearOp.hpp>
71 
72 //#include <Thyra_LinearOpBase.hpp>
73 #ifdef HAVE_XPETRA_TPETRA
74 //#include <Thyra_TpetraLinearOp.hpp>
75 #include <Thyra_TpetraThyraWrappers.hpp>
76 #include <Thyra_TpetraVector.hpp>
77 #include <Tpetra_Map.hpp>
78 #include <Tpetra_Vector.hpp>
79 #include <Tpetra_CrsMatrix.hpp>
80 #include <Xpetra_TpetraMap.hpp>
81 #include <Xpetra_TpetraCrsMatrix.hpp>
82 #endif
83 #ifdef HAVE_XPETRA_EPETRA
84 #include <Thyra_EpetraLinearOp.hpp>
85 #include <Thyra_EpetraThyraWrappers.hpp>
86 #include <Thyra_SpmdVectorBase.hpp>
87 #include <Thyra_get_Epetra_Operator.hpp>
88 #include <Epetra_Map.h>
89 #include <Epetra_Vector.h>
90 #include <Epetra_CrsMatrix.h>
91 #include <Xpetra_EpetraMap.hpp>
92 #include <Xpetra_EpetraCrsMatrix.hpp>
93 #endif
94 
95 namespace Xpetra {
96 
97 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> class BlockedCrsMatrix;
98 
99 template <class Scalar,
100  class LocalOrdinal = int,
101  class GlobalOrdinal = LocalOrdinal,
102  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
103 class ThyraUtils {
104 
105 private:
106 #undef XPETRA_THYRAUTILS_SHORT
107 #include "Xpetra_UseShortNames.hpp"
108 
109 public:
110 
111  static Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >
112  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) {
113 
114  Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map = toXpetra(vectorSpace);
115 
116  if(stridedBlockId == -1) {
117  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo.size() != 0);
118  }
119  else {
120  TEUCHOS_TEST_FOR_EXCEPT(map->getNodeNumElements() % stridingInfo[stridedBlockId] != 0);
121  }
122 
123  Teuchos::RCP<Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> > ret = Xpetra::StridedMapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(map, stridingInfo, stridedBlockId, offset);
124  return ret;
125  }
126 
127  static Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
128  toXpetra(const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >& vectorSpace, const Teuchos::RCP<const Teuchos::Comm<int> >& comm) {
129 
130  // check whether we have a Tpetra based Thyra operator
131  bool bIsTpetra = false;
132  #ifdef HAVE_XPETRA_TPETRA
133  Teuchos::RCP<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_vsc = Teuchos::rcp_dynamic_cast<const Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vectorSpace);
134  bIsTpetra = Teuchos::is_null(tpetra_vsc) ? false : true;
135  #endif
136 
137  // check whether we have an Epetra based Thyra operator
138  bool bIsEpetra = !bIsTpetra; // note: this is a little bit fragile!
139 
140  #ifdef HAVE_XPETRA_TPETRA
141  if(bIsTpetra) {
142  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
143  //Teuchos::RCP<const Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rgVec = vectorSpace->createMember();
144  Teuchos::RCP<Thyra::VectorBase<Scalar> > rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string("label"));
145  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgVec));
146  Teuchos::RCP<const Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rgTpetraVec = TOE::getTpetraVector(rgVec);
147  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraVec));
148  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rgTpetraMap = rgTpetraVec->getMap();
149  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgTpetraMap));
150 
151  Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node > > rgXpetraMap = Xpetra::toXpetraNonConst(rgTpetraMap);
152  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgXpetraMap));
153  return rgXpetraMap;
154  }
155  #endif
156 
157  #ifdef HAVE_XPETRA_EPETRA
158  if(bIsEpetra) {
159  /*const Teuchos::RCP<const Thyra::SpmdVectorSpaceBase<double> > spmd_vs =
160  Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double> >(vectorSpace, true);
161  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(spmd_vs));*/
162  //Teuchos::RCP<const Epetra_Comm> epComm = Thyra::get_Epetra_Comm(*spmd_vs->getComm());
163  //Teuchos::RCP<const Epetra_Comm> epComm = Thyra::get_Epetra_Comm(*comm);
164  Teuchos::RCP<const Epetra_Comm> epComm = Xpetra::toEpetra(comm);
165 
166  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epComm));
167  //Teuchos::RCP<const Epetra_Comm> epComm = Thyra::get_Epetra_Comm(const Teuchos::Comm<Ordinal>& comm);
168  Teuchos::RCP<const Epetra_Map> rgEpetraMap = Thyra::get_Epetra_Map(*vectorSpace, epComm);
169  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgEpetraMap));
170 
171  Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node > > rgXpetraMap = Teuchos::rcp(new Xpetra::EpetraMap(rgEpetraMap));
172  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(rgXpetraMap));
173  return rgXpetraMap;
174  }
175  #endif
176  return Teuchos::null;
177  }
178 
179  static bool isTpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
180  // check whether we have a Tpetra based Thyra operator
181  bool bIsTpetra = false;
182 #ifdef HAVE_XPETRA_TPETRA
183  Teuchos::RCP<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tpetra_op = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(op);
184  bIsTpetra = Teuchos::is_null(tpetra_op) ? false : true;
185 #endif
186  return bIsTpetra;
187  }
188 
189  static bool isEpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > & op){
190  // check whether we have an Epetra based Thyra operator
191  bool bIsEpetra = false;
192 #ifdef HAVE_XPETRA_EPETRA
193  Teuchos::RCP<const Thyra::EpetraLinearOp> epetra_op = Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op);
194  bIsEpetra = Teuchos::is_null(epetra_op) ? false : true;
195 #endif
196  return bIsEpetra;
197  }
198 
199  static Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
200  toXpetra(const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >& op) {
201 
202 #ifdef HAVE_XPETRA_TPETRA
203  if(isTpetra(op)) {
204  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
205  Teuchos::RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getConstTpetraOperator(op);
206  // we should also add support for the const versions!
207  //getConstTpetraOperator(const RCP<const LinearOpBase<Scalar> > &op);
208  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
209  Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
210  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraRowMat));
211  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
212  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraCrsMat));
213  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraNcnstCrsMat = Teuchos::rcp_const_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraCrsMat);
214  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraNcnstCrsMat));
215 
216  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xTpetraCrsMat =
217  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(TpetraNcnstCrsMat));
218  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
219 
220  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ret =
221  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xTpetraCrsMat);
222  return ret;
223  }
224 #endif
225 
226 #ifdef HAVE_XPETRA_EPETRA
227  if(isEpetra(op)) {
228  Teuchos::RCP<const Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator( *op );
229  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
230  Teuchos::RCP<const Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_op);
231  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_rowmat));
232  Teuchos::RCP<const Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(epetra_rowmat);
233  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_crsmat));
234  Teuchos::RCP<Epetra_CrsMatrix> epetra_ncnstcrsmat = Teuchos::rcp_const_cast<Epetra_CrsMatrix>(epetra_crsmat);
235  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_ncnstcrsmat));
236 
237  Teuchos::RCP<Xpetra::EpetraCrsMatrix > xEpetraCrsMat =
238  Teuchos::rcp(new Xpetra::EpetraCrsMatrix(epetra_ncnstcrsmat));
239  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
240 
241  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > ret =
242  Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xEpetraCrsMat);
243  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ret));
244  return ret;
245  }
246 #endif
247  return Teuchos::null;
248  }
249 
250  static Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
251  toXpetra(const Teuchos::RCP<Thyra::LinearOpBase<Scalar> >& op) {
252 
253 #ifdef HAVE_XPETRA_TPETRA
254  if(isTpetra(op)) {
255  typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
256  Teuchos::RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraOp = TOE::getTpetraOperator(op);
257  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraOp));
258  Teuchos::RCP<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraOp);
259  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraRowMat));
260  Teuchos::RCP<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > TpetraCrsMat = Teuchos::rcp_dynamic_cast<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(TpetraRowMat);
261  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(TpetraCrsMat));
262 
263  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xTpetraCrsMat =
264  Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(TpetraCrsMat));
265  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpetraCrsMat));
266  return xTpetraCrsMat;
267  }
268 #endif
269 
270 #ifdef HAVE_XPETRA_EPETRA
271  if(isEpetra(op)) {
272  Teuchos::RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator( *op );
273  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_op));
274  Teuchos::RCP<Epetra_RowMatrix> epetra_rowmat = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetra_op);
275  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_rowmat));
276  Teuchos::RCP<Epetra_CrsMatrix> epetra_crsmat = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetra_rowmat);
277  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epetra_crsmat));
278 
279  Teuchos::RCP<Xpetra::EpetraCrsMatrix > xEpetraCrsMat =
280  Teuchos::rcp(new Xpetra::EpetraCrsMatrix(epetra_crsmat));
281  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpetraCrsMat));
282  return xEpetraCrsMat;
283  }
284 #endif
285  return Teuchos::null;
286  }
287 
288  static Teuchos::RCP<const Thyra::LinearOpBase<Scalar> >
289  toThyra(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& mat) {
290  // create a Thyra operator from Xpetra::CrsMatrix
291  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
292 
293 #ifdef HAVE_XPETRA_TPETRA
294  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
295  if(tpetraMat!=Teuchos::null) {
296  Teuchos::RCP<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xTpCrsMat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
297  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpCrsMat));
298  Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrix();
299  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
300 
301  Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
302  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpRowMat));
303  Teuchos::RCP<const Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
304  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpOperator));
305 
306  thyraOp = Thyra::createConstLinearOp(tpOperator);
307  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
308  }
309 #endif
310 
311 #ifdef HAVE_XPETRA_EPETRA
312  Teuchos::RCP<const Xpetra::EpetraCrsMatrix> epetraMat = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraCrsMatrix>(mat);
313  if(epetraMat!=Teuchos::null) {
314  Teuchos::RCP<Xpetra::EpetraCrsMatrix > xEpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrix >(mat);
315  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpCrsMat));
316  Teuchos::RCP<const Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrix();
317  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
318 
319  Teuchos::RCP<const Thyra::EpetraLinearOp> thyraEpOp = Thyra::epetraLinearOp(epCrsMat,"op");
320  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
321  thyraOp = thyraEpOp;
322  }
323 #endif
324  return thyraOp;
325  }
326 
327  static Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
328  toThyra(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& mat) {
329  // create a Thyra operator from Xpetra::CrsMatrix
330  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thyraOp = Teuchos::null;
331 
332 #ifdef HAVE_XPETRA_TPETRA
333  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
334  if(tpetraMat!=Teuchos::null) {
335  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xTpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
336  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xTpCrsMat));
337  Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpCrsMat = xTpCrsMat->getTpetra_CrsMatrixNonConst();
338  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpCrsMat));
339 
340  Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpRowMat = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpCrsMat);
341  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpRowMat));
342  Teuchos::RCP<Tpetra::Operator <Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOperator = Teuchos::rcp_dynamic_cast<Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(tpRowMat);
343  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpOperator));
344 
345  thyraOp = Thyra::createLinearOp(tpOperator);
346  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraOp));
347  }
348 #endif
349 
350 #ifdef HAVE_XPETRA_EPETRA
351  Teuchos::RCP<Xpetra::EpetraCrsMatrix> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrix>(mat);
352  if(epetraMat!=Teuchos::null) {
353  Teuchos::RCP<Xpetra::EpetraCrsMatrix > xEpCrsMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrix >(mat);
354  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xEpCrsMat));
355  Teuchos::RCP<Epetra_CrsMatrix> epCrsMat = xEpCrsMat->getEpetra_CrsMatrixNonConst();
356  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(epCrsMat));
357 
358  Teuchos::RCP<Thyra::EpetraLinearOp> thyraEpOp = Thyra::nonconstEpetraLinearOp(epCrsMat,"op");
359  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraEpOp));
360  thyraOp = thyraEpOp;
361  }
362 #endif
363  return thyraOp;
364  }
365 
366  static Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
367  toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& mat) {
368 
369  int nRows = mat->Rows();
370  int nCols = mat->Cols();
371 
372  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ablock = mat->getMatrix(0,0);
373 
374  bool bTpetra = false;
375  bool bEpetra = false;
376 #ifdef HAVE_XPETRA_TPETRA
377  Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(Ablock);
378  if(tpetraMat!=Teuchos::null) bTpetra = true;
379 #endif
380 
381 #ifdef HAVE_XPETRA_EPETRA
382  Teuchos::RCP<Xpetra::EpetraCrsMatrix> epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrix>(Ablock);
383  if(epetraMat!=Teuchos::null) bEpetra = true;
384  #endif
385 
386  TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra
387 
388  // create new Thyra blocked operator
389  Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> > blockMat =
390  Thyra::defaultBlockedLinearOp<Scalar>();
391 
392  blockMat->beginBlockFill(nRows,nCols);
393 
394  for (int r=0; r<nRows; ++r) {
395  for (int c=0; c<nCols; ++c) {
396  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thBlock =
397  Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(mat->getMatrix(r,c));
398  std::stringstream label; label << "A" << r << c;
399  blockMat->setBlock(r,c,thBlock);
400  }
401  }
402 
403  blockMat->endBlockFill();
404 
405  return blockMat;
406  }
407 
408 }; // end Utils class
409 
410 } // end namespace Xpetra
411 
412 #define XPETRA_THYRAUTILS_SHORT
413 #endif // HAVE_XPETRA_THYRA
414 
415 #endif // XPETRA_THYRAUTILS_HPP
Xpetra::toEpetra
const Epetra_CrsGraph & toEpetra(const RCP< const CrsGraph< int, GlobalOrdinal > > &graph)
Definition: Xpetra_EpetraCrsGraph.cpp:57
Xpetra_TpetraMap.hpp
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
Xpetra_UseShortNames.hpp
Xpetra
Xpetra namespace
Definition: Xpetra_CrsGraph.hpp:63
Xpetra::toXpetra
RCP< const CrsGraph< int, GlobalOrdinal > > toXpetra(const Epetra_CrsGraph &g)
Definition: Xpetra_EpetraCrsGraph.cpp:159
Xpetra_StridedMapFactory.hpp
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra::StridedMapFactory::Build
static RCP< StridedMap > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
Definition: Xpetra_StridedMapFactory.hpp:88
Xpetra_StridedMap.hpp
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrixNonConst
RCP< Epetra_CrsMatrix > getEpetra_CrsMatrixNonConst() const
Get the underlying Epetra matrix.
Definition: Xpetra_EpetraCrsMatrix.hpp:369
Xpetra::toXpetraNonConst
const RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > toXpetraNonConst(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &)
Definition: Xpetra_TpetraMap.hpp:297
Xpetra_EpetraCrsMatrix.hpp
Xpetra_EpetraMap.hpp
Xpetra::EpetraCrsMatrixT::getEpetra_CrsMatrix
RCP< const Epetra_CrsMatrix > getEpetra_CrsMatrix() const
Get the underlying Epetra matrix.
Definition: Xpetra_EpetraCrsMatrix.hpp:366
Xpetra_ConfigDefs.hpp
Xpetra_MapExtractor.hpp
Xpetra_Map.hpp
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Xpetra_TpetraCrsMatrix.hpp
Generated by   doxygen 1.8.11