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