Skip to content

Commit

Permalink
Fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@cgcgcg
cgcgcg committed May 23, 2020
1 parent 6e164ad commit 1848f5f
Show file tree
Hide file tree
Showing 2 changed files with 58 additions and 78 deletions.
22 changes: 0 additions & 22 deletions packages/muelu/adapters/xpetra/MueLu_RefMaxwell_decl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -410,28 +410,6 @@ namespace MueLu {
xpExporter->setDistributorParameters(matvecParams);
}

//! find nonzeros in ArrayRCP
void FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
Teuchos::ArrayRCP<bool> nonzeros);

//!
void DetectDirichletCols(const Matrix& A,
const Teuchos::ArrayRCP<const bool>& dirichletRows,
Teuchos::ArrayRCP<bool> dirichletCols,
Teuchos::ArrayRCP<bool> dirichletDomain);

#ifdef HAVE_MUELU_KOKKOS_REFACTOR
//! find nonzeros in Kokkos::View
void FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_um vals,
Kokkos::View<bool*, typename Node::device_type> nonzeros);

//!
void DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Kokkos::View<const bool*, typename Node::device_type> & dirichletRows,
Kokkos::View<bool*, typename Node::device_type> dirichletCols,
Kokkos::View<bool*, typename Node::device_type> dirichletDomain);
#endif

//! Two hierarchies: one for the coarse (1,1)-block, another for the (2,2)-block
Teuchos::RCP<Hierarchy> HierarchyH_, Hierarchy22_;
Teuchos::RCP<SmootherBase> PreSmoother_, PostSmoother_;
Expand Down
114 changes: 58 additions & 56 deletions packages/muelu/adapters/xpetra/MueLu_RefMaxwell_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,56 +106,30 @@
namespace MueLu {

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getDomainMap() const {
return SM_Matrix_->getDomainMap();
}


template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getRangeMap() const {
return SM_Matrix_->getRangeMap();
}

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
Teuchos::ArrayRCP<bool> nonzeros) {
void FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
Teuchos::ArrayRCP<bool> nonzeros) {
TEUCHOS_ASSERT(vals.size() == nonzeros.size());
typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
const magnitudeType eps = 2.0*Teuchos::ScalarTraits<magnitudeType>::eps();
for(size_t i=0; i<static_cast<size_t>(vals.size()); i++) {
nonzeros[i] = (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > eps);
}
}

#ifdef HAVE_MUELU_KOKKOS_REFACTOR
template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_um vals,
Kokkos::View<bool*, typename Node::device_type> nonzeros) {
using ATS = Kokkos::ArithTraits<Scalar>;
using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
TEUCHOS_ASSERT(vals.extent(0) == nonzeros.extent(0));
const typename ATS::magnitudeType eps = 2.0*ATS::eps();

Kokkos::parallel_for("MueLu:RefMaxwell::FindNonZeros", range_type(0,vals.extent(0)),
KOKKOS_LAMBDA (const size_t i) {
nonzeros(i) = (ATS::magnitude(vals(i,0)) > eps);
});
}
#endif


template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Teuchos::ArrayRCP<const bool>& dirichletRows,
Teuchos::ArrayRCP<bool> dirichletCols,
Teuchos::ArrayRCP<bool> dirichletDomain) {
void DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Teuchos::ArrayRCP<const bool>& dirichletRows,
Teuchos::ArrayRCP<bool> dirichletCols,
Teuchos::ArrayRCP<bool> dirichletDomain) {
const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
RCP<const Map> domMap = A.getDomainMap();
RCP<const Map> rowMap = A.getRowMap();
RCP<const Map> colMap = A.getColMap();
RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == rowMap->getNodeNumElements());
TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == colMap->getNodeNumElements());
TEUCHOS_ASSERT(static_cast<size_t>(dirichletDomain.size()) == domMap->getNodeNumElements());
RCP<MultiVector> myColsToZero = MultiVectorFactory::Build(colMap, 1, /*zeroOut=*/true);
RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, 1, /*zeroOut=*/true);
// Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
if (dirichletRows[i]) {
Expand All @@ -167,10 +141,10 @@ namespace MueLu {
}
}

RCP<MultiVector> globalColsToZero;
RCP<const Import> importer = A.getCrsGraph()->getImporter();
RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
if (!importer.is_null()) {
globalColsToZero = MultiVectorFactory::Build(domMap, 1, /*zeroOut=*/true);
globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, 1, /*zeroOut=*/true);
// export to domain map
globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
// import to column map
Expand All @@ -179,30 +153,46 @@ namespace MueLu {
else
globalColsToZero = myColsToZero;

FindNonZeros(globalColsToZero->getData(0),dirichletDomain);
FindNonZeros(myColsToZero->getData(0),dirichletCols);
FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(globalColsToZero->getData(0),dirichletDomain);
FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(myColsToZero->getData(0),dirichletCols);
}


#ifdef HAVE_MUELU_KOKKOS_REFACTOR

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_um vals,
Kokkos::View<bool*, typename Node::device_type> nonzeros) {
using ATS = Kokkos::ArithTraits<Scalar>;
using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
TEUCHOS_ASSERT(vals.extent(0) == nonzeros.extent(0));
const typename ATS::magnitudeType eps = 2.0*ATS::eps();

Kokkos::parallel_for("MueLu:RefMaxwell::FindNonZeros", range_type(0,vals.extent(0)),
KOKKOS_LAMBDA (const size_t i) {
nonzeros(i) = (ATS::magnitude(vals(i,0)) > eps);
});
}

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Kokkos::View<const bool*, typename Node::device_type> & dirichletRows,
Kokkos::View<bool*, typename Node::device_type> dirichletCols,
Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
void DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Kokkos::View<const bool*, typename Node::device_type> & dirichletRows,
Kokkos::View<bool*, typename Node::device_type> dirichletCols,
Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
RCP<const Map> domMap = A.getDomainMap();
RCP<const Map> rowMap = A.getRowMap();
RCP<const Map> colMap = A.getColMap();
RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
TEUCHOS_ASSERT(dirichletRows.extent(0) == rowMap->getNodeNumElements());
TEUCHOS_ASSERT(dirichletCols.extent(0) == colMap->getNodeNumElements());
TEUCHOS_ASSERT(dirichletDomain.extent(0) == domMap->getNodeNumElements());
RCP<Vector> myColsToZero = VectorFactory::Build(colMap, /*zeroOut=*/true);
RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, /*zeroOut=*/true);
// Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
auto myColsToZeroView = myColsToZero->template getLocalView<typename Node::device_type>();
auto localMatrix = A.getLocalMatrix();
Kokkos::parallel_for("MueLu:RefMaxwell::DetectDirichletCols", range_type(0,rowMap->getNodeNumElements()),
KOKKOS_LAMBDA(const LO row) {
KOKKOS_LAMBDA(const LocalOrdinal row) {
if (dirichletRows(row)) {
auto rowView = localMatrix.row(row);
auto length = rowView.length;
Expand All @@ -212,23 +202,35 @@ namespace MueLu {
}
});

RCP<Vector> globalColsToZero;
RCP<const Import> importer = A.getCrsGraph()->getImporter();
RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
if (!importer.is_null()) {
globalColsToZero = VectorFactory::Build(domMap, /*zeroOut=*/true);
globalColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, /*zeroOut=*/true);
// export to domain map
globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
// import to column map
myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
}
else
globalColsToZero = myColsToZero;
FindNonZeros(globalColsToZero->template getLocalView<typename Node::device_type>(),dirichletDomain);
FindNonZeros(myColsToZero->template getLocalView<typename Node::device_type>(),dirichletCols);
FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(globalColsToZero->template getLocalView<typename Node::device_type>(),dirichletDomain);
FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(myColsToZero->template getLocalView<typename Node::device_type>(),dirichletCols);
}

#endif

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getDomainMap() const {
return SM_Matrix_->getDomainMap();
}


template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getRangeMap() const {
return SM_Matrix_->getRangeMap();
}


template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameters(Teuchos::ParameterList& list) {

Expand Down

0 comments on commit 1848f5f

Please sign in to comment.