47 #ifndef XPETRA_THYRAUTILS_HPP
48 #define XPETRA_THYRAUTILS_HPP
51 #ifdef HAVE_XPETRA_THYRA
55 #ifdef HAVE_XPETRA_TPETRA
56 #include "Tpetra_ConfigDefs.hpp"
59 #ifdef HAVE_XPETRA_EPETRA
60 #include "Epetra_config.h"
61 #include "Epetra_CombineMode.h"
73 #include <Thyra_VectorSpaceBase.hpp>
74 #include <Thyra_SpmdVectorSpaceBase.hpp>
75 #include <Thyra_ProductVectorSpaceBase.hpp>
76 #include <Thyra_ProductMultiVectorBase.hpp>
77 #include <Thyra_VectorSpaceBase.hpp>
78 #include <Thyra_DefaultProductVectorSpace.hpp>
79 #include <Thyra_DefaultBlockedLinearOp.hpp>
80 #include <Thyra_LinearOpBase.hpp>
81 #include <Thyra_DetachedMultiVectorView.hpp>
82 #include <Thyra_MultiVectorStdOps.hpp>
84 #ifdef HAVE_XPETRA_TPETRA
85 #include <Thyra_TpetraThyraWrappers.hpp>
86 #include <Thyra_TpetraVector.hpp>
87 #include <Thyra_TpetraVectorSpace.hpp>
88 #include <Tpetra_Map.hpp>
89 #include <Tpetra_Vector.hpp>
90 #include <Tpetra_CrsMatrix.hpp>
92 #include <Xpetra_TpetraCrsMatrix.hpp>
94 #ifdef HAVE_XPETRA_EPETRA
95 #include <Thyra_EpetraLinearOp.hpp>
96 #include <Thyra_EpetraThyraWrappers.hpp>
97 #include <Thyra_SpmdVectorBase.hpp>
98 #include <Thyra_get_Epetra_Operator.hpp>
99 #include <Epetra_Map.h>
100 #include <Epetra_Vector.h>
101 #include <Epetra_CrsMatrix.h>
108 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
class BlockedCrsMatrix;
110 template <
class Scalar,
111 class LocalOrdinal = int,
112 class GlobalOrdinal = LocalOrdinal,
117 #undef XPETRA_THYRAUTILS_SHORT
126 if(stridedBlockId == -1) {
140 using Teuchos::rcp_dynamic_cast;
142 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
143 typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
145 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
147 RCP<Map> resultMap = Teuchos::null;
148 RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<
const ThyProdVecSpaceBase >(vectorSpace);
149 if(prodVectorSpace != Teuchos::null) {
152 TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error,
"Found a product vector space with zero blocks.");
153 std::vector<RCP<const Map> > mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
154 std::vector<RCP<const Map> > mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
155 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
156 RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
163 std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
164 for(
int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
165 gidOffsets[i] = mapsThyra[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
168 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
169 RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
176 #ifdef HAVE_XPETRA_TPETRA
182 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
183 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
184 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
185 typedef Thyra::VectorBase<Scalar> ThyVecBase;
186 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string(
"label"));
188 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
190 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
194 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.");
205 using Teuchos::rcp_dynamic_cast;
207 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
208 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
212 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
215 RCP<MultiVector> xpMultVec = Teuchos::null;
219 if(thyProdVec != Teuchos::null) {
228 RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec);
232 for (
int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
233 RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
236 xpBlockedMultVec->setMultiVector(b, xpBlockMV,
true );
240 #ifdef HAVE_XPETRA_TPETRA
241 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
242 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
244 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
245 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
247 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<
const ThySpmdMultVecBase>(v);
248 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<
const ThyTpMultVec>(thyraSpmdMultiVector);
250 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
252 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
254 xpMultVec =
rcp(
new XpTpMultVec(tpNonConstMultVec));
256 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.");
267 Teuchos::rcp_const_cast<
const Thyra::MultiVectorBase<Scalar> >(v);
274 static bool isTpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
278 bool bIsTpetra =
false;
279 #ifdef HAVE_XPETRA_TPETRA
285 #ifdef HAVE_XPETRA_EPETRA
286 Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
288 Teuchos::rcp_dynamic_cast<
const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
290 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
291 if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
292 std::cout <<
"ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
293 std::cout <<
" If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
294 std::cout <<
" properly set!" << std::endl;
295 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op,
true) << std::endl;
304 if(bIsTpetra ==
false) {
306 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
307 if(ThyBlockedOp != Teuchos::null) {
310 ThyBlockedOp->getBlock(0,0);
312 bIsTpetra = isTpetra(b00);
320 static bool isEpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
324 static bool isBlockedOperator(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
327 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
328 if(ThyBlockedOp != Teuchos::null) {
337 #ifdef HAVE_XPETRA_TPETRA
339 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
361 #ifdef HAVE_XPETRA_EPETRA
366 return Teuchos::null;
372 #ifdef HAVE_XPETRA_TPETRA
374 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
385 return xTpetraCrsMat;
389 #ifdef HAVE_XPETRA_EPETRA
394 return Teuchos::null;
402 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(map);
403 if(bmap.is_null() ==
false) {
406 for(
size_t i = 0; i < bmap->getNumMaps(); i++) {
409 Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(bmap->getMap(i,
true));
413 thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
418 #ifdef HAVE_XPETRA_TPETRA
421 if (tpetraMap == Teuchos::null)
422 throw Exceptions::BadCast(
"Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
423 RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetraMap->getTpetra_Map());
424 thyraMap = thyraTpetraMap;
428 #ifdef HAVE_XPETRA_EPETRA
442 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getGlobalNumElements())!=thMap->dim(), std::logic_error,
"Global dimension of Xpetra map and Thyra VectorSpaceBase are different.");
444 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error,
"Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
445 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getNodeNumElements())!=thSpmdMap->localSubDim(), std::logic_error,
"Local dimension of Xpetra map and Thyra VectorSpaceBase on one (or more) processor(s) are different.");
450 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMVec == Teuchos::null, std::logic_error,
"Cannot cast MultiVectorBase to SpmdMultiVectorBase.");
453 const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
454 const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
459 for(
size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
462 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
463 (*thyData)(i,j) = vecData[i];
475 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getGlobalNumElements())!=thMap->dim(), std::logic_error,
"Global dimension of Xpetra map and Thyra VectorSpaceBase are different.");
477 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error,
"Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
478 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getNodeNumElements())!=thSpmdMap->localSubDim(), std::logic_error,
"Local dimension of Xpetra map and Thyra VectorSpaceBase on one (or more) processor(s) are different.");
486 const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
487 const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
492 for(
size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
495 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
496 (*thyData)(i,j) = vecData[i];
508 using Teuchos::rcp_dynamic_cast;
510 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
511 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
512 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
515 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
519 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
520 if(prodTarget != Teuchos::null) {
521 RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<
const BlockedMultiVector>(source);
522 if(bSourceVec.is_null() ==
true) {
526 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error,
"Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
527 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.");
529 for(
int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
531 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb,
false);
535 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
536 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(vs);
537 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
538 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
539 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
543 for(
size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
547 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
548 (*thyData)(i,j) = xpData[i];
555 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.");
557 for(
int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
559 RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb,
true);
565 Thyra::assign(thySubBlock.
ptr(), *thyXpSubBlock);
574 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(target->range());
575 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error,
"Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
576 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
577 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
578 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
582 for(
size_t j = 0; j < source->getNumVectors(); ++j) {
585 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
586 (*thyData)(i,j) = xpData[i];
599 #ifdef HAVE_XPETRA_TPETRA
601 if(tpetraMat!=Teuchos::null) {
613 thyraOp = Thyra::createConstLinearOp(tpOperator);
616 #ifdef HAVE_XPETRA_EPETRA
636 #ifdef HAVE_XPETRA_TPETRA
638 if(tpetraMat!=Teuchos::null) {
650 thyraOp = Thyra::createLinearOp(tpOperator);
654 #ifdef HAVE_XPETRA_EPETRA
670 int nRows = mat->Rows();
671 int nCols = mat->Cols();
677 #ifdef HAVE_XPETRA_TPETRA
679 if(tpetraMat!=Teuchos::null) {
683 Thyra::defaultBlockedLinearOp<Scalar>();
685 blockMat->beginBlockFill(nRows,nCols);
687 for (
int r=0; r<nRows; ++r) {
688 for (
int c=0; c<nCols; ++c) {
691 if(xpmat == Teuchos::null)
continue;
698 if(xpblock != Teuchos::null) {
699 if(xpblock->Rows() == 1 && xpblock->Cols() == 1) {
703 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpwrap->getCrsMatrix());
706 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpblock);
712 thBlock = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpwrap->getCrsMatrix());
715 blockMat->setBlock(r,c,thBlock);
719 blockMat->endBlockFill();
724 #ifdef HAVE_XPETRA_EPETRA
731 #endif // endif HAVE_XPETRA_TPETRA
735 #if defined(HAVE_XPETRA_EPETRA) && !defined(HAVE_XPETRA_TPETRA)
738 #endif // endif HAVE_XPETRA_EPETRA
746 #ifdef HAVE_XPETRA_EPETRA
748 #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
750 class ThyraUtils<double, int, int,
EpetraNode> {
752 typedef double Scalar;
753 typedef int LocalOrdinal;
754 typedef int GlobalOrdinal;
758 #undef XPETRA_THYRAUTILS_SHORT
768 if(stridedBlockId == -1) {
782 using Teuchos::rcp_dynamic_cast;
784 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
785 typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
788 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
790 RCP<Map> resultMap = Teuchos::null;
792 RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<
const ThyProdVecSpaceBase >(vectorSpace);
793 if(prodVectorSpace != Teuchos::null) {
796 TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error,
"Found a product vector space with zero blocks.");
797 std::vector<RCP<const Map> > mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
798 std::vector<RCP<const Map> > mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
799 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
800 RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
807 std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
808 for(
int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
809 gidOffsets[i] = mapsThyra[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
812 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
813 RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
824 bool bIsTpetra =
false;
825 #ifdef HAVE_XPETRA_TPETRA
826 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
827 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
834 bool bIsEpetra = !bIsTpetra;
836 #ifdef HAVE_XPETRA_TPETRA
838 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
839 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
840 typedef Thyra::VectorBase<Scalar> ThyVecBase;
841 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
842 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
843 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
844 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string(
"label"));
846 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
848 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
858 #ifdef HAVE_XPETRA_EPETRA
861 RCP<const Epetra_Map>
862 epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map> >(vectorSpace,
"epetra_map");
880 using Teuchos::rcp_dynamic_cast;
882 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
883 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
887 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
890 RCP<MultiVector> xpMultVec = Teuchos::null;
894 if(thyProdVec != Teuchos::null) {
903 RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec);
907 for (
int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
908 RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
911 xpBlockedMultVec->setMultiVector(b, xpBlockMV,
true );
919 bool bIsTpetra =
false;
920 #ifdef HAVE_XPETRA_TPETRA
921 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
922 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
926 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
927 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
928 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
930 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
932 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<
const ThySpmdMultVecBase>(v);
933 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<
const ThyTpMultVec>(thyraSpmdMultiVector);
934 if(thyraTpetraMultiVector != Teuchos::null) {
936 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
938 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
940 xpMultVec =
rcp(
new XpTpMultVec(tpNonConstMultVec));
945 #ifdef HAVE_XPETRA_EPETRA
946 if(bIsTpetra ==
false) {
952 RCP<const Epetra_Map> eMap = Teuchos::rcpFromRef(xeMap->getEpetra_Map());
956 RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<
Epetra_MultiVector>(epMultVec);
971 Teuchos::rcp_const_cast<
const Thyra::MultiVectorBase<Scalar> >(v);
977 static bool isTpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
979 bool bIsTpetra =
false;
980 #ifdef HAVE_XPETRA_TPETRA
981 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
982 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
989 #ifdef HAVE_XPETRA_EPETRA
990 Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
992 Teuchos::rcp_dynamic_cast<
const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
994 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
995 if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
996 std::cout <<
"ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
997 std::cout <<
" If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
998 std::cout <<
" properly set!" << std::endl;
999 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op,
true) << std::endl;
1009 if(bIsTpetra ==
false) {
1011 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
1012 if(ThyBlockedOp != Teuchos::null) {
1015 ThyBlockedOp->getBlock(0,0);
1017 bIsTpetra = isTpetra(b00);
1025 static bool isEpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
1027 bool bIsEpetra =
false;
1029 #ifdef HAVE_XPETRA_EPETRA
1038 if(bIsEpetra ==
false) {
1040 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op,
false);
1041 if(ThyBlockedOp != Teuchos::null) {
1044 ThyBlockedOp->getBlock(0,0);
1046 bIsEpetra = isEpetra(b00);
1054 static bool isBlockedOperator(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
1057 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
1058 if(ThyBlockedOp != Teuchos::null) {
1067 #ifdef HAVE_XPETRA_TPETRA
1069 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1070 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1072 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1097 #ifdef HAVE_XPETRA_EPETRA
1118 return Teuchos::null;
1124 #ifdef HAVE_XPETRA_TPETRA
1126 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1127 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1129 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1140 return xTpetraCrsMat;
1147 #ifdef HAVE_XPETRA_EPETRA
1159 return xEpetraCrsMat;
1162 return Teuchos::null;
1170 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(map);
1171 if(bmap.is_null() ==
false) {
1174 for(
size_t i = 0; i < bmap->getNumMaps(); i++) {
1177 Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(bmap->getMap(i,
true));
1181 thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
1186 #ifdef HAVE_XPETRA_TPETRA
1188 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1189 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1191 if (tpetraMap == Teuchos::null)
1192 throw Exceptions::BadCast(
"Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
1193 RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > tpMap = tpetraMap->getTpetra_Map();
1194 RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
1195 thyraMap = thyraTpetraMap;
1202 #ifdef HAVE_XPETRA_EPETRA
1205 if (epetraMap == Teuchos::null)
1206 throw Exceptions::BadCast(
"Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
1207 RCP<const Thyra::VectorSpaceBase<Scalar> > thyraEpetraMap = Thyra::create_VectorSpace(Teuchos::rcpFromRef(epetraMap->getEpetra_Map()));
1208 thyraMap = thyraEpetraMap;
1220 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getGlobalNumElements())!=thMap->dim(), std::logic_error,
"Global dimension of Xpetra map and Thyra VectorSpaceBase are different.");
1222 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error,
"Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
1223 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getNodeNumElements())!=thSpmdMap->localSubDim(), std::logic_error,
"Local dimension of Xpetra map and Thyra VectorSpaceBase on one (or more) processor(s) are different.");
1228 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMVec == Teuchos::null, std::logic_error,
"Cannot cast MultiVectorBase to SpmdMultiVectorBase.");
1231 const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
1232 const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
1237 for(
size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
1240 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1241 (*thyData)(i,j) = vecData[i];
1253 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getGlobalNumElements())!=thMap->dim(), std::logic_error,
"Global dimension of Xpetra map and Thyra VectorSpaceBase are different.");
1255 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error,
"Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
1256 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getNodeNumElements())!=thSpmdMap->localSubDim(), std::logic_error,
"Local dimension of Xpetra map and Thyra VectorSpaceBase on one (or more) processor(s) are different.");
1264 const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
1265 const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
1270 for(
size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
1273 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1274 (*thyData)(i,j) = vecData[i];
1283 using Teuchos::rcp_dynamic_cast;
1285 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1286 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
1287 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1290 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1294 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
1295 if(prodTarget != Teuchos::null) {
1297 RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<
const BlockedMultiVector>(source);
1298 if(bSourceVec.is_null() ==
true) {
1302 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error,
"Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
1303 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.");
1305 for(
int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1307 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb,
false);
1311 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
1312 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(vs);
1313 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
1314 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
1315 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1319 for(
size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
1323 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1324 (*thyData)(i,j) = xpData[i];
1331 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.");
1333 for(
int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1335 RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb,
true);
1341 Thyra::assign(thySubBlock.
ptr(), *thyXpSubBlock);
1350 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(target->range());
1351 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error,
"Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
1352 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
1353 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
1354 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1358 for(
size_t j = 0; j < source->getNumVectors(); ++j) {
1361 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1362 (*thyData)(i,j) = xpData[i];
1373 #ifdef HAVE_XPETRA_TPETRA
1375 if(tpetraMat!=Teuchos::null) {
1376 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1377 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1389 thyraOp = Thyra::createConstLinearOp(tpOperator);
1397 #ifdef HAVE_XPETRA_EPETRA
1399 if(epetraMat!=Teuchos::null) {
1407 thyraOp = thyraEpOp;
1418 #ifdef HAVE_XPETRA_TPETRA
1420 if(tpetraMat!=Teuchos::null) {
1421 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1422 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
1434 thyraOp = Thyra::createLinearOp(tpOperator);
1442 #ifdef HAVE_XPETRA_EPETRA
1444 if(epetraMat!=Teuchos::null) {
1452 thyraOp = thyraEpOp;
1462 #endif // #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
1464 #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
1466 class ThyraUtils<double, int, long long,
EpetraNode> {
1468 typedef double Scalar;
1469 typedef int LocalOrdinal;
1470 typedef long long GlobalOrdinal;
1474 #undef XPETRA_THYRAUTILS_SHORT
1484 if(stridedBlockId == -1) {
1498 using Teuchos::rcp_dynamic_cast;
1500 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1501 typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
1503 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
1505 RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<
const ThyProdVecSpaceBase >(vectorSpace);
1506 if(prodVectorSpace != Teuchos::null) {
1509 TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error,
"Found a product vector space with zero blocks.");
1510 std::vector<RCP<const Map> > mapsThyra(prodVectorSpace->numBlocks(), Teuchos::null);
1511 std::vector<RCP<const Map> > mapsXpetra(prodVectorSpace->numBlocks(), Teuchos::null);
1512 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
1513 RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
1520 std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
1521 for(
int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
1522 gidOffsets[i] = mapsThyra[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
1525 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
1526 RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
1539 bool bIsTpetra =
false;
1540 #ifdef HAVE_XPETRA_TPETRA
1541 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1542 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1549 bool bIsEpetra = !bIsTpetra;
1551 #ifdef HAVE_XPETRA_TPETRA
1553 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1554 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1555 typedef Thyra::VectorBase<Scalar> ThyVecBase;
1556 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
1557 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
1558 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1559 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string(
"label"));
1561 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
1563 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
1575 #ifdef HAVE_XPETRA_EPETRA
1578 RCP<const Epetra_Map>
1579 epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map> >(vectorSpace,
"epetra_map");
1598 using Teuchos::rcp_dynamic_cast;
1600 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1601 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1605 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
1608 RCP<MultiVector> xpMultVec = Teuchos::null;
1612 if(thyProdVec != Teuchos::null) {
1621 RCP<BlockedMultiVector> xpBlockedMultVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xpMultVec);
1625 for (
int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
1626 RCP<const ThyMultVecBase> thyBlockMV = thyProdVec->getMultiVectorBlock(b);
1629 xpBlockedMultVec->setMultiVector(b, xpBlockMV,
true );
1637 bool bIsTpetra =
false;
1638 #ifdef HAVE_XPETRA_TPETRA
1639 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1640 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1644 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
1645 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
1646 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
1648 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
1650 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<
const ThySpmdMultVecBase>(v);
1651 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<
const ThyTpMultVec>(thyraSpmdMultiVector);
1652 if(thyraTpetraMultiVector != Teuchos::null) {
1654 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
1656 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
1658 xpMultVec =
rcp(
new XpTpMultVec(tpNonConstMultVec));
1665 #ifdef HAVE_XPETRA_EPETRA
1666 if(bIsTpetra ==
false) {
1672 RCP<const Epetra_Map> eMap = Teuchos::rcpFromRef(xeMap->getEpetra_Map());
1676 RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<
Epetra_MultiVector>(epMultVec);
1692 Teuchos::rcp_const_cast<
const Thyra::MultiVectorBase<Scalar> >(v);
1698 static bool isTpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
1700 bool bIsTpetra =
false;
1701 #ifdef HAVE_XPETRA_TPETRA
1702 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1703 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1710 #ifdef HAVE_XPETRA_EPETRA
1711 Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
1713 Teuchos::rcp_dynamic_cast<
const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
1715 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
1716 if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
1717 std::cout <<
"ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
1718 std::cout <<
" If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
1719 std::cout <<
" properly set!" << std::endl;
1720 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op,
true) << std::endl;
1730 if(bIsTpetra ==
false) {
1732 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
1733 if(ThyBlockedOp != Teuchos::null) {
1736 ThyBlockedOp->getBlock(0,0);
1738 bIsTpetra = isTpetra(b00);
1746 static bool isEpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
1748 bool bIsEpetra =
false;
1750 #ifdef HAVE_XPETRA_EPETRA
1759 if(bIsEpetra ==
false) {
1761 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op,
false);
1762 if(ThyBlockedOp != Teuchos::null) {
1765 ThyBlockedOp->getBlock(0,0);
1767 bIsEpetra = isEpetra(b00);
1775 static bool isBlockedOperator(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
1778 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
1779 if(ThyBlockedOp != Teuchos::null) {
1788 #ifdef HAVE_XPETRA_TPETRA
1790 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1791 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1793 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1818 #ifdef HAVE_XPETRA_EPETRA
1839 return Teuchos::null;
1845 #ifdef HAVE_XPETRA_TPETRA
1847 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1848 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1850 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1861 return xTpetraCrsMat;
1868 #ifdef HAVE_XPETRA_EPETRA
1880 return xEpetraCrsMat;
1883 return Teuchos::null;
1891 RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(map);
1892 if(bmap.is_null() ==
false) {
1895 for(
size_t i = 0; i < bmap->getNumMaps(); i++) {
1898 Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(bmap->getMap(i,
true));
1902 thyraMap = Thyra::productVectorSpace<Scalar>(vecSpaces());
1907 #ifdef HAVE_XPETRA_TPETRA
1909 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
1910 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
1912 if (tpetraMap == Teuchos::null)
1913 throw Exceptions::BadCast(
"Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
1914 RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > tpMap = tpetraMap->getTpetra_Map();
1915 RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
1916 thyraMap = thyraTpetraMap;
1923 #ifdef HAVE_XPETRA_EPETRA
1926 if (epetraMap == Teuchos::null)
1927 throw Exceptions::BadCast(
"Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
1928 RCP<const Thyra::VectorSpaceBase<Scalar> > thyraEpetraMap = Thyra::create_VectorSpace(Teuchos::rcpFromRef(epetraMap->getEpetra_Map()));
1929 thyraMap = thyraEpetraMap;
1941 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getGlobalNumElements())!=thMap->dim(), std::logic_error,
"Global dimension of Xpetra map and Thyra VectorSpaceBase are different.");
1943 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error,
"Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
1944 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getNodeNumElements())!=thSpmdMap->localSubDim(), std::logic_error,
"Local dimension of Xpetra map and Thyra VectorSpaceBase on one (or more) processor(s) are different.");
1949 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMVec == Teuchos::null, std::logic_error,
"Cannot cast MultiVectorBase to SpmdMultiVectorBase.");
1952 const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
1953 const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
1958 for(
size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
1961 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1962 (*thyData)(i,j) = vecData[i];
1974 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getGlobalNumElements())!=thMap->dim(), std::logic_error,
"Global dimension of Xpetra map and Thyra VectorSpaceBase are different.");
1976 TEUCHOS_TEST_FOR_EXCEPTION(thSpmdMap == Teuchos::null, std::logic_error,
"Cannot cast VectorSpaceBase to SpmdVectorSpaceBase.");
1977 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Teuchos::Ordinal>(vec->getMap()->getNodeNumElements())!=thSpmdMap->localSubDim(), std::logic_error,
"Local dimension of Xpetra map and Thyra VectorSpaceBase on one (or more) processor(s) are different.");
1985 const LocalOrdinal localOffset = ( thSpmdMap != Teuchos::null ? thSpmdMap->localOffset() : 0 );
1986 const LocalOrdinal localSubDim = ( thSpmdMap != Teuchos::null ? thSpmdMap->localSubDim() : thMap->dim() );
1991 for(
size_t j = 0; j < Teuchos::as<size_t>(thSpmdMVec->domain()->dim()); ++j) {
1994 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1995 (*thyData)(i,j) = vecData[i];
2004 using Teuchos::rcp_dynamic_cast;
2006 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
2007 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
2008 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
2011 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
2015 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
2016 if(prodTarget != Teuchos::null) {
2017 RCP<const BlockedMultiVector> bSourceVec = rcp_dynamic_cast<
const BlockedMultiVector>(source);
2018 if(bSourceVec.is_null() ==
true) {
2022 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error,
"Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
2023 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.");
2025 for(
int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
2027 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb,
false);
2031 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
2032 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(vs);
2033 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
2034 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
2035 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
2039 for(
size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
2043 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
2044 (*thyData)(i,j) = xpData[i];
2051 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.");
2053 for(
int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
2055 RCP<MultiVector> xpSubBlock = bSourceVec->getMultiVector(bbb,
true);
2061 Thyra::assign(thySubBlock.
ptr(), *thyXpSubBlock);
2070 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(target->range());
2071 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error,
"Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
2072 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
2073 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
2074 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
2078 for(
size_t j = 0; j < source->getNumVectors(); ++j) {
2081 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
2082 (*thyData)(i,j) = xpData[i];
2093 #ifdef HAVE_XPETRA_TPETRA
2095 if(tpetraMat!=Teuchos::null) {
2096 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
2097 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
2109 thyraOp = Thyra::createConstLinearOp(tpOperator);
2117 #ifdef HAVE_XPETRA_EPETRA
2119 if(epetraMat!=Teuchos::null) {
2127 thyraOp = thyraEpOp;
2138 #ifdef HAVE_XPETRA_TPETRA
2140 if(tpetraMat!=Teuchos::null) {
2141 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
2142 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_LONG_LONG) && defined(HAVE_TPETRA_INST_DOUBLE)))
2154 thyraOp = Thyra::createLinearOp(tpOperator);
2162 #ifdef HAVE_XPETRA_EPETRA
2164 if(epetraMat!=Teuchos::null) {
2172 thyraOp = thyraEpOp;
2182 #endif // XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
2184 #endif // HAVE_XPETRA_EPETRA
2188 #define XPETRA_THYRAUTILS_SHORT
2189 #endif // HAVE_XPETRA_THYRA
2191 #endif // XPETRA_THYRAUTILS_HPP
static RCP< StridedMap > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
Map constructor with Xpetra-defined contiguous uniform distribution.
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Exception throws to report errors in the internal logical of the program.
Xpetra utility class for common map-related routines.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
const RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > toXpetraNonConst(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &)
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
TypeTo as(const TypeFrom &t)
Concrete implementation of Xpetra::Matrix.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)