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