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" 69 #include <Thyra_VectorSpaceBase.hpp> 70 #include <Thyra_ProductVectorSpaceBase.hpp> 71 #include <Thyra_ProductMultiVectorBase.hpp> 72 #include <Thyra_VectorSpaceBase.hpp> 73 #include <Thyra_DefaultBlockedLinearOp.hpp> 74 #include <Thyra_LinearOpBase.hpp> 75 #include <Thyra_DetachedMultiVectorView.hpp> 77 #ifdef HAVE_XPETRA_TPETRA 78 #include <Thyra_TpetraThyraWrappers.hpp> 79 #include <Thyra_TpetraVector.hpp> 80 #include <Thyra_TpetraVectorSpace.hpp> 81 #include <Tpetra_Map.hpp> 82 #include <Tpetra_Vector.hpp> 83 #include <Tpetra_CrsMatrix.hpp> 87 #ifdef HAVE_XPETRA_EPETRA 88 #include <Thyra_EpetraLinearOp.hpp> 89 #include <Thyra_EpetraThyraWrappers.hpp> 90 #include <Thyra_SpmdVectorBase.hpp> 91 #include <Thyra_get_Epetra_Operator.hpp> 92 #include <Epetra_Map.h> 93 #include <Epetra_Vector.h> 94 #include <Epetra_CrsMatrix.h> 101 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
class BlockedCrsMatrix;
103 template <
class Scalar,
104 class LocalOrdinal = int,
105 class GlobalOrdinal = LocalOrdinal,
110 #undef XPETRA_THYRAUTILS_SHORT 120 if(stridedBlockId == -1) {
134 using Teuchos::rcp_dynamic_cast;
136 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
142 typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
145 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
147 RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<
const ThyProdVecSpaceBase >(vectorSpace);
148 if(prodVectorSpace != Teuchos::null) {
151 TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error,
"Found a product vector space with zero blocks.");
152 std::vector<RCP<Map> > maps(prodVectorSpace->numBlocks(), Teuchos::null);
153 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
154 RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
160 std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
161 for(
int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
162 gidOffsets[i] = maps[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
166 std::vector<GlobalOrdinal> gids;
167 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
168 for(LocalOrdinal l = 0; l < as<LocalOrdinal>(maps[b]->getNodeNumElements()); ++l) {
169 GlobalOrdinal gid = maps[b]->getGlobalElement(l) + gidOffsets[b];
177 RCP<Map> fullMap =
MapFactory::Build(maps[0]->lib(), INVALID, gidsView, maps[0]->getIndexBase(), comm);
182 #ifdef HAVE_XPETRA_TPETRA 188 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
189 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
190 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
191 typedef Thyra::VectorBase<Scalar> ThyVecBase;
192 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string(
"label"));
194 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
196 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
203 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.");
207 return Teuchos::null;
214 using Teuchos::rcp_dynamic_cast;
216 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
217 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
225 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
230 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
235 RCP<MultiVector> xpMultVec = Teuchos::null;
239 if(thyProdVec != Teuchos::null) {
248 std::vector<GlobalOrdinal> lidOffsets(thyProdVec->productSpace()->numBlocks()+1,0);
249 for (
int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
252 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
253 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : thySubMap->dim() );
254 lidOffsets[b+1] = localSubDim + lidOffsets[b];
255 RCP<Thyra::ConstDetachedMultiVectorView<Scalar> > thyData =
256 Teuchos::rcp(
new Thyra::ConstDetachedMultiVectorView<Scalar>(thyProdVec->getMultiVectorBlock(b),
Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
257 for(
size_t vv = 0; vv < xpMultVec->getNumVectors(); ++vv) {
258 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
259 xpMultVec->replaceLocalValue(i + lidOffsets[b] , vv, (*thyData)(i,vv));
267 #ifdef HAVE_XPETRA_TPETRA 268 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
269 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
271 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
272 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
274 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<
const ThySpmdMultVecBase>(v);
275 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<
const ThyTpMultVec>(thyraSpmdMultiVector);
277 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
279 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
281 xpMultVec =
rcp(
new XpTpMultVec(tpNonConstMultVec));
285 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.");
289 return Teuchos::null;
296 Teuchos::rcp_const_cast<
const Thyra::MultiVectorBase<Scalar> >(v);
303 static bool isTpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
307 bool bIsTpetra =
false;
308 #ifdef HAVE_XPETRA_TPETRA 314 #ifdef HAVE_XPETRA_EPETRA
315 Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
317 Teuchos::rcp_dynamic_cast<
const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
319 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
320 if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
321 std::cout <<
"ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
322 std::cout <<
" If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
323 std::cout <<
" properly set!" << std::endl;
324 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op,
true) << std::endl;
333 if(bIsTpetra ==
false) {
335 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
336 if(ThyBlockedOp != Teuchos::null) {
339 ThyBlockedOp->getBlock(0,0);
341 bIsTpetra = isTpetra(b00);
349 static bool isEpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
353 static bool isBlockedOperator(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
356 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
357 if(ThyBlockedOp != Teuchos::null) {
366 #ifdef HAVE_XPETRA_TPETRA 368 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
390 #ifdef HAVE_XPETRA_EPETRA 395 return Teuchos::null;
401 #ifdef HAVE_XPETRA_TPETRA 403 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
414 return xTpetraCrsMat;
418 #ifdef HAVE_XPETRA_EPETRA 423 return Teuchos::null;
429 #ifdef HAVE_XPETRA_TPETRA 432 if (tpetraMap == Teuchos::null)
433 throw Exceptions::BadCast(
"Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
434 RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetraMap->getTpetra_Map());
435 thyraMap = thyraTpetraMap;
439 #ifdef HAVE_XPETRA_EPETRA 453 using Teuchos::rcp_dynamic_cast;
455 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
456 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
457 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
460 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
464 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
465 if(prodTarget != Teuchos::null) {
469 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error,
"Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
470 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.");
472 for(
int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
474 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb,
false);
478 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
479 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(vs);
480 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
481 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
482 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
486 for(
size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
490 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
491 (*thyData)(i,j) = xpData[i];
500 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(target->range());
501 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error,
"Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
502 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
503 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
504 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
508 for(
size_t j = 0; j < source->getNumVectors(); ++j) {
511 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
512 (*thyData)(i,j) = xpData[i];
523 bool bIsTpetra =
false;
525 #ifdef HAVE_XPETRA_TPETRA 527 if(tpetraMat!=Teuchos::null) {
539 thyraOp = Thyra::createConstLinearOp(tpOperator);
544 #ifdef HAVE_XPETRA_EPETRA 555 bool bIsTpetra =
false;
557 #ifdef HAVE_XPETRA_TPETRA 559 if(tpetraMat!=Teuchos::null) {
571 thyraOp = Thyra::createLinearOp(tpOperator);
576 #ifdef HAVE_XPETRA_EPETRA 585 int nRows = mat->Rows();
586 int nCols = mat->Cols();
590 bool bTpetra =
false;
591 #ifdef HAVE_XPETRA_TPETRA 593 if(tpetraMat!=Teuchos::null) bTpetra =
true;
596 #ifdef HAVE_XPETRA_EPETRA 602 Thyra::defaultBlockedLinearOp<Scalar>();
604 blockMat->beginBlockFill(nRows,nCols);
606 for (
int r=0; r<nRows; ++r) {
607 for (
int c=0; c<nCols; ++c) {
609 Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(mat->getMatrix(r,c));
610 std::stringstream label; label <<
"A" << r << c;
611 blockMat->setBlock(r,c,thBlock);
615 blockMat->endBlockFill();
625 #ifdef HAVE_XPETRA_EPETRA 627 class ThyraUtils<double, int, int,
EpetraNode> {
629 typedef double Scalar;
630 typedef int LocalOrdinal;
631 typedef int GlobalOrdinal;
635 #undef XPETRA_THYRAUTILS_SHORT 645 if(stridedBlockId == -1) {
659 using Teuchos::rcp_dynamic_cast;
661 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
662 typedef Thyra::ProductVectorSpaceBase<Scalar> ThyProdVecSpaceBase;
665 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
667 RCP<const ThyProdVecSpaceBase > prodVectorSpace = rcp_dynamic_cast<
const ThyProdVecSpaceBase >(vectorSpace);
668 if(prodVectorSpace != Teuchos::null) {
671 TEUCHOS_TEST_FOR_EXCEPTION(prodVectorSpace->numBlocks()==0, std::logic_error,
"Found a product vector space with zero blocks.");
672 std::vector<RCP<Map> > maps(prodVectorSpace->numBlocks(), Teuchos::null);
673 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
674 RCP<const ThyVecSpaceBase > bv = prodVectorSpace->getBlock(b);
680 std::vector<GlobalOrdinal> gidOffsets(prodVectorSpace->numBlocks(),0);
681 for(
int i = 1; i < prodVectorSpace->numBlocks(); ++i) {
682 gidOffsets[i] = maps[i-1]->getMaxAllGlobalIndex() + gidOffsets[i-1] + 1;
686 std::vector<GlobalOrdinal> gids;
687 for (
int b = 0; b < prodVectorSpace->numBlocks(); ++b){
688 for(LocalOrdinal l = 0; l < as<LocalOrdinal>(maps[b]->getNodeNumElements()); ++l) {
689 GlobalOrdinal gid = maps[b]->getGlobalElement(l) + gidOffsets[b];
697 RCP<Map> fullMap =
MapFactory::Build(maps[0]->lib(), INVALID, gidsView, maps[0]->getIndexBase(), comm);
706 bool bIsTpetra =
false;
707 #ifdef HAVE_XPETRA_TPETRA 708 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 709 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 716 bool bIsEpetra = !bIsTpetra;
718 #ifdef HAVE_XPETRA_TPETRA 720 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 721 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 722 typedef Thyra::VectorBase<Scalar> ThyVecBase;
723 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpMap;
724 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpVector;
725 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
726 RCP<ThyVecBase> rgVec = Thyra::createMember<Scalar>(vectorSpace, std::string(
"label"));
728 RCP<const TpVector> rgTpetraVec = TOE::getTpetraVector(rgVec);
730 RCP<const TpMap> rgTpetraMap = rgTpetraVec->getMap();
742 #ifdef HAVE_XPETRA_EPETRA 745 RCP<const Epetra_Map> epMap = Teuchos::null;
746 RCP<const Epetra_Map>
747 epetra_map = Teuchos::get_extra_data<RCP<const Epetra_Map> >(vectorSpace,
"epetra_map");
759 return Teuchos::null;
766 using Teuchos::rcp_dynamic_cast;
768 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
769 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
772 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
777 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyUtils;
780 RCP<MultiVector> xpMultVec = Teuchos::null;
784 if(thyProdVec != Teuchos::null) {
793 std::vector<GlobalOrdinal> lidOffsets(thyProdVec->productSpace()->numBlocks()+1,0);
794 for (
int b = 0; b < thyProdVec->productSpace()->numBlocks(); ++b){
797 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
798 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : thySubMap->dim() );
799 lidOffsets[b+1] = localSubDim + lidOffsets[b];
800 RCP<Thyra::ConstDetachedMultiVectorView<Scalar> > thyData =
801 Teuchos::rcp(
new Thyra::ConstDetachedMultiVectorView<Scalar>(thyProdVec->getMultiVectorBlock(b),
Teuchos::Range1D(localOffset,localOffset+localSubDim-1)));
802 for(
size_t vv = 0; vv < xpMultVec->getNumVectors(); ++vv) {
803 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
804 xpMultVec->replaceLocalValue(i + lidOffsets[b] , vv, (*thyData)(i,vv));
813 bool bIsTpetra =
false;
814 #ifdef HAVE_XPETRA_TPETRA 815 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 816 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 820 typedef Thyra::SpmdMultiVectorBase<Scalar> ThySpmdMultVecBase;
821 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> ConverterT;
822 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpMultVec;
824 typedef Thyra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyTpMultVec;
826 RCP<const ThySpmdMultVecBase> thyraSpmdMultiVector = rcp_dynamic_cast<
const ThySpmdMultVecBase>(v);
827 RCP<const ThyTpMultVec> thyraTpetraMultiVector = rcp_dynamic_cast<
const ThyTpMultVec>(thyraSpmdMultiVector);
828 if(thyraTpetraMultiVector != Teuchos::null) {
830 const RCP<const TpMultVec> tpMultVec = ConverterT::getConstTpetraMultiVector(v);
832 RCP<TpMultVec> tpNonConstMultVec = Teuchos::rcp_const_cast<TpMultVec>(tpMultVec);
834 xpMultVec =
rcp(
new XpTpMultVec(tpNonConstMultVec));
841 #ifdef HAVE_XPETRA_EPETRA 842 if(bIsTpetra ==
false) {
848 RCP<const Epetra_Map> eMap = Teuchos::rcpFromRef(xeMap->getEpetra_Map());
852 RCP<Epetra_MultiVector> epNonConstMultVec = Teuchos::rcp_const_cast<Epetra_MultiVector>(epMultVec);
861 return Teuchos::null;
868 Teuchos::rcp_const_cast<
const Thyra::MultiVectorBase<Scalar> >(v);
874 static bool isTpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
876 bool bIsTpetra =
false;
877 #ifdef HAVE_XPETRA_TPETRA 878 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 879 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 886 #ifdef HAVE_XPETRA_EPETRA
887 Teuchos::rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op) == Teuchos::null &&
889 Teuchos::rcp_dynamic_cast<
const Thyra::DefaultBlockedLinearOp<Scalar> >(op) == Teuchos::null) {
891 typedef Thyra::TpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraLinearOp_t;
892 if(Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op) == Teuchos::null) {
893 std::cout <<
"ATTENTION: The dynamic cast to the TpetraLinearOp failed even though it should be a TpetraLinearOp." << std::endl;
894 std::cout <<
" If you are using Panzer or Stratimikos you might check that the template parameters are " << std::endl;
895 std::cout <<
" properly set!" << std::endl;
896 std::cout << Teuchos::rcp_dynamic_cast<const TpetraLinearOp_t>(op,
true) << std::endl;
906 if(bIsTpetra ==
false) {
908 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
909 if(ThyBlockedOp != Teuchos::null) {
912 ThyBlockedOp->getBlock(0,0);
914 bIsTpetra = isTpetra(b00);
922 static bool isEpetra(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
924 bool bIsEpetra =
false;
926 #ifdef HAVE_XPETRA_EPETRA 935 if(bIsEpetra ==
false) {
937 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op,
false);
938 if(ThyBlockedOp != Teuchos::null) {
941 ThyBlockedOp->getBlock(0,0);
943 bIsEpetra = isEpetra(b00);
951 static bool isBlockedOperator(
const Teuchos::RCP<
const Thyra::LinearOpBase<Scalar> > & op){
954 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<Scalar> >(op);
955 if(ThyBlockedOp != Teuchos::null) {
964 #ifdef HAVE_XPETRA_TPETRA 966 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 967 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 969 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
994 #ifdef HAVE_XPETRA_EPETRA 1015 return Teuchos::null;
1021 #ifdef HAVE_XPETRA_TPETRA 1023 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 1024 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 1026 typedef Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node> TOE;
1037 return xTpetraCrsMat;
1044 #ifdef HAVE_XPETRA_EPETRA 1056 return xEpetraCrsMat;
1059 return Teuchos::null;
1065 #ifdef HAVE_XPETRA_TPETRA 1067 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 1068 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 1070 if (tpetraMap == Teuchos::null)
1071 throw Exceptions::BadCast(
"Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::TpetraMap failed");
1072 RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > tpMap = tpetraMap->getTpetra_Map();
1073 RCP<Thyra::TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thyraTpetraMap = Thyra::tpetraVectorSpace<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpMap);
1074 thyraMap = thyraTpetraMap;
1081 #ifdef HAVE_XPETRA_EPETRA 1084 if (epetraMap == Teuchos::null)
1085 throw Exceptions::BadCast(
"Xpetra::ThyraUtils::toThyra: Cast from Xpetra::Map to Xpetra::EpetraMap failed");
1086 RCP<const Thyra::VectorSpaceBase<Scalar> > thyraEpetraMap = Thyra::create_VectorSpace(Teuchos::rcpFromRef(epetraMap->getEpetra_Map()));
1087 thyraMap = thyraEpetraMap;
1096 using Teuchos::rcp_dynamic_cast;
1098 typedef Thyra::VectorSpaceBase<Scalar> ThyVecSpaceBase;
1099 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThySpmdVecSpaceBase;
1100 typedef Thyra::MultiVectorBase<Scalar> ThyMultVecBase;
1103 typedef Thyra::ProductMultiVectorBase<Scalar> ThyProdMultVecBase;
1107 RCP<ThyProdMultVecBase> prodTarget = rcp_dynamic_cast<ThyProdMultVecBase>(target);
1108 if(prodTarget != Teuchos::null) {
1112 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor == Teuchos::null, std::logic_error,
"Found a Thyra product vector, but user did not provide an Xpetra::MapExtractor.");
1113 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.");
1115 for(
int bbb = 0; bbb < prodTarget->productSpace()->numBlocks(); ++bbb) {
1117 RCP<MultiVector> xpSubBlock = mapExtractor->ExtractVector(source, bbb,
false);
1121 RCP<const ThyVecSpaceBase> vs = thySubBlock->range();
1122 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(vs);
1123 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
1124 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : vs->dim() );
1125 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1129 for(
size_t j = 0; j < xpSubBlock->getNumVectors(); ++j) {
1133 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1134 (*thyData)(i,j) = xpData[i];
1143 RCP<const ThySpmdVecSpaceBase> mpi_vs = rcp_dynamic_cast<
const ThySpmdVecSpaceBase>(target->range());
1144 TEUCHOS_TEST_FOR_EXCEPTION(mpi_vs == Teuchos::null, std::logic_error,
"Failed to cast Thyra::VectorSpaceBase to Thyra::SpmdVectorSpaceBase.");
1145 const LocalOrdinal localOffset = ( mpi_vs != Teuchos::null ? mpi_vs->localOffset() : 0 );
1146 const LocalOrdinal localSubDim = ( mpi_vs != Teuchos::null ? mpi_vs->localSubDim() : target->range()->dim() );
1147 RCP<Thyra::DetachedMultiVectorView<Scalar> > thyData =
1151 for(
size_t j = 0; j < source->getNumVectors(); ++j) {
1154 for(LocalOrdinal i = 0; i < localSubDim; ++i) {
1155 (*thyData)(i,j) = xpData[i];
1166 #ifdef HAVE_XPETRA_TPETRA 1168 if(tpetraMat!=Teuchos::null) {
1169 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 1170 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 1182 thyraOp = Thyra::createConstLinearOp(tpOperator);
1190 #ifdef HAVE_XPETRA_EPETRA 1192 if(epetraMat!=Teuchos::null) {
1200 thyraOp = thyraEpOp;
1211 #ifdef HAVE_XPETRA_TPETRA 1213 if(tpetraMat!=Teuchos::null) {
1214 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 1215 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 1227 thyraOp = Thyra::createLinearOp(tpOperator);
1235 #ifdef HAVE_XPETRA_EPETRA 1237 if(epetraMat!=Teuchos::null) {
1245 thyraOp = thyraEpOp;
1256 int nRows = mat->Rows();
1257 int nCols = mat->Cols();
1261 bool bTpetra =
false;
1262 bool bEpetra =
false;
1263 #ifdef HAVE_XPETRA_TPETRA 1264 #if ((defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \ 1265 (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE))) 1267 if(tpetraMat!=Teuchos::null) bTpetra =
true;
1273 #ifdef HAVE_XPETRA_EPETRA 1275 if(epetraMat!=Teuchos::null) bEpetra =
true;
1282 Thyra::defaultBlockedLinearOp<Scalar>();
1284 blockMat->beginBlockFill(nRows,nCols);
1286 for (
int r=0; r<nRows; ++r) {
1287 for (
int c=0; c<nCols; ++c) {
1289 Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(mat->getMatrix(r,c));
1290 std::stringstream label; label <<
"A" << r << c;
1291 blockMat->setBlock(r,c,thBlock);
1295 blockMat->endBlockFill();
1302 #endif // HAVE_XPETRA_EPETRA 1305 #define XPETRA_THYRAUTILS_SHORT 1306 #endif // HAVE_XPETRA_THYRA 1308 #endif // XPETRA_THYRAUTILS_HPP static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
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.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static RCP< StridedMap > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
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)
TypeTo as(const TypeFrom &t)
Create an Xpetra::Map instance.
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)