Skip to content

Commit

Permalink
Merge pull request #7227 from cgcgcg/muelubugs
Browse files Browse the repository at this point in the history
Xpetra: Fix for CheckRepairMainDiagonal
  • Loading branch information
cgcgcg authored Apr 22, 2020
2 parents 485357a + 5045938 commit f86a556
Show file tree
Hide file tree
Showing 3 changed files with 43 additions and 8 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -1353,14 +1353,17 @@ namespace MueLu {
TEUCHOS_TEST_FOR_EXCEPTION(partName != "zoltan" && partName != "zoltan2", Exceptions::InvalidArgument,
"Invalid partitioner name: \"" << partName << "\". Valid options: \"zoltan\", \"zoltan2\"");

bool switched = false;
(void)switched;
#ifndef HAVE_MUELU_ZOLTAN
bool switched = false;
if (partName == "zoltan") {
this->GetOStream(Warnings0) << "Zoltan interface is not available, trying to switch to Zoltan2" << std::endl;
partName = "zoltan2";
switched = true;
}
#else
# ifndef HAVE_MUELU_ZOLTAN
bool switched = false;
# endif
#endif
#ifndef HAVE_MUELU_ZOLTAN2
if (partName == "zoltan2" && !switched) {
Expand Down
40 changes: 36 additions & 4 deletions packages/muelu/test/scaling/DriverCore.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -223,6 +223,40 @@ void PreconditionerSetup(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalO
}


#if defined(HAVE_MUELU_EPETRA)

// Helper functions for compilation purposes
template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
struct Matvec_Wrapper{
static void UnwrapEpetra(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& X,
Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& B,
Teuchos::RCP<const Epetra_CrsMatrix>& Aepetra,
Teuchos::RCP<Epetra_MultiVector>& Xepetra,
Teuchos::RCP<Epetra_MultiVector>& Bepetra) {
throw std::runtime_error("Template parameter mismatch");
}
};


template<class GlobalOrdinal>
struct Matvec_Wrapper<double,int,GlobalOrdinal,Kokkos::Compat::KokkosSerialWrapperNode> {
static void UnwrapEpetra(Teuchos::RCP<Xpetra::Matrix<double,int,GlobalOrdinal,Kokkos::Compat::KokkosSerialWrapperNode> >& A,
Teuchos::RCP<Xpetra::MultiVector<double,int,GlobalOrdinal,Kokkos::Compat::KokkosSerialWrapperNode> >& X,
Teuchos::RCP<Xpetra::MultiVector<double,int,GlobalOrdinal,Kokkos::Compat::KokkosSerialWrapperNode> >& B,
Teuchos::RCP<const Epetra_CrsMatrix>& Aepetra,
Teuchos::RCP<Epetra_MultiVector>& Xepetra,
Teuchos::RCP<Epetra_MultiVector>& Bepetra) {
typedef double SC;
typedef int LO;
typedef GlobalOrdinal GO;
typedef Kokkos::Compat::KokkosSerialWrapperNode NO;
Aepetra = Xpetra::Helpers<SC, LO, GO, NO>::Op2EpetraCrs(A);
Xepetra = Teuchos::rcp(& Xpetra::toEpetra(*X),false);
Bepetra = Teuchos::rcp(& Xpetra::toEpetra(*B),false);
}
};
#endif

//*************************************************************************************
template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Expand Down Expand Up @@ -271,13 +305,11 @@ void SystemSolve(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,N
Btpetra = rcp(& Xpetra::toTpetra(*B),false);
}
#endif
#if defined(HAVE_MUELU_EPETRA) && !defined(HAVE_MUELU_INST_COMPLEX_INT_INT) && !defined(HAVE_MUELU_INST_FLOAT_INT_INT)
#if defined(HAVE_MUELU_EPETRA)
Teuchos::RCP<const Epetra_CrsMatrix> Aepetra;
Teuchos::RCP<Epetra_MultiVector> Xepetra,Bepetra;
if(lib==Xpetra::UseEpetra) {
Aepetra = Xpetra::Helpers<SC, LO, GO, NO>::Op2EpetraCrs(A);
Xepetra = rcp(& Xpetra::toEpetra(*X),false);
Bepetra = rcp(& Xpetra::toEpetra(*B),false);
Matvec_Wrapper<SC,LO,GO,NO>::UnwrapEpetra(A,X,B,Aepetra,Xepetra,Bepetra);
}
#endif

Expand Down
4 changes: 2 additions & 2 deletions packages/xpetra/sup/Utils/Xpetra_MatrixUtils.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -453,7 +453,7 @@ class MatrixUtils {
Ac->getLocalDiagCopy(*diagVec);

LocalOrdinal lZeroDiags = 0;
Teuchos::ArrayRCP< Scalar > diagVal = diagVec->getDataNonConst(0);
Teuchos::ArrayRCP< const Scalar > diagVal = diagVec->getData(0);

for (size_t i = 0; i < rowMap->getNodeNumElements(); i++) {
if (TST::magnitude(diagVal[i]) <= threshold) {
Expand Down Expand Up @@ -521,7 +521,7 @@ class MatrixUtils {
#ifdef HAVE_XPETRA_DEBUG // only for debugging
// check whether Ac has been repaired...
Ac->getLocalDiagCopy(*diagVec);
Teuchos::ArrayRCP< Scalar > diagVal2 = diagVec->getDataNonConst(0);
diagVal = diagVec->getData(0);
for (size_t r = 0; r < Ac->getRowMap()->getNodeNumElements(); r++) {
if (TST::magnitude(diagVal[r]) <= threshold) {
fos << "Error: there are too small entries left on diagonal after repair..." << std::endl;
Expand Down

0 comments on commit f86a556

Please sign in to comment.