46 #ifndef MUELU_UTILITIES_KOKKOS_DEF_HPP 47 #define MUELU_UTILITIES_KOKKOS_DEF_HPP 50 #include <Teuchos_DefaultComm.hpp> 51 #include <Teuchos_ParameterList.hpp> 55 #ifdef HAVE_MUELU_EPETRA 57 # include "Epetra_MpiComm.h" 61 #include <Kokkos_Core.hpp> 62 #include <KokkosSparse_CrsMatrix.hpp> 63 #include <KokkosSparse_getDiagCopy.hpp> 65 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT) 72 #include <Xpetra_EpetraUtils.hpp> 73 #include <Xpetra_EpetraMultiVector.hpp> 77 #ifdef HAVE_MUELU_TPETRA 78 #include <MatrixMarket_Tpetra.hpp> 79 #include <Tpetra_RowMatrixTransposer.hpp> 80 #include <TpetraExt_MatrixMatrix.hpp> 81 #include <Xpetra_TpetraMultiVector.hpp> 82 #include <Xpetra_TpetraCrsMatrix.hpp> 83 #include <Xpetra_TpetraBlockCrsMatrix.hpp> 86 #ifdef HAVE_MUELU_EPETRA 87 #include <Xpetra_EpetraMap.hpp> 90 #include <Xpetra_BlockedCrsMatrix.hpp> 91 #include <Xpetra_DefaultPlatform.hpp> 92 #include <Xpetra_Import.hpp> 93 #include <Xpetra_ImportFactory.hpp> 94 #include <Xpetra_Map.hpp> 95 #include <Xpetra_MapFactory.hpp> 96 #include <Xpetra_Matrix.hpp> 97 #include <Xpetra_MatrixMatrix.hpp> 98 #include <Xpetra_MatrixFactory.hpp> 99 #include <Xpetra_MultiVector.hpp> 100 #include <Xpetra_MultiVectorFactory.hpp> 101 #include <Xpetra_Operator.hpp> 102 #include <Xpetra_Vector.hpp> 103 #include <Xpetra_VectorFactory.hpp> 107 #include <KokkosKernels_Handle.hpp> 108 #include <KokkosGraph_RCM.hpp> 114 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
115 Teuchos::RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
116 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
117 GetMatrixDiagonal(
const Matrix& A) {
118 const auto rowMap = A.getRowMap();
119 auto diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap,
true);
121 A.getLocalDiagCopy(*diag);
126 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
127 Teuchos::RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
128 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
129 GetMatrixDiagonalInverse(
const Matrix& A, Magnitude tol,
const bool doLumped) {
130 Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer(
"Utilities_kokkos::GetMatrixDiagonalInverse");
132 using local_matrix_type =
typename Matrix::local_matrix_type;
133 using local_graph_type =
typename local_matrix_type::staticcrsgraph_type;
134 using value_type =
typename local_matrix_type::value_type;
135 using ordinal_type =
typename local_matrix_type::ordinal_type;
136 using execution_space =
typename local_matrix_type::execution_space;
137 using memory_space =
typename local_matrix_type::memory_space;
143 using KAT = Kokkos::ArithTraits<value_type>;
146 RCP<const Map> rowMap = A.getRowMap();
147 RCP<Vector> diag = VectorFactory::Build(rowMap,
false);
150 local_matrix_type localMatrix = A.getLocalMatrixDevice();
151 auto diagVals = diag->getDeviceLocalView(Xpetra::Access::ReadWrite);
153 ordinal_type numRows = localMatrix.graph.numRows();
160 Kokkos::parallel_for(
"Utilities_kokkos::GetMatrixDiagonalInverse",
162 KOKKOS_LAMBDA(
const ordinal_type rowIdx) {
163 bool foundDiagEntry =
false;
164 auto myRow = localMatrix.rowConst(rowIdx);
165 for(ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
166 if(myRow.colidx(entryIdx) == rowIdx) {
167 foundDiagEntry =
true;
168 if(KAT::magnitude(myRow.value(entryIdx)) > KAT::magnitude(tol)) {
169 diagVals(rowIdx, 0) = KAT::one() / myRow.value(entryIdx);
171 diagVals(rowIdx, 0) = KAT::zero();
177 if(!foundDiagEntry) {diagVals(rowIdx, 0) = KAT::zero();}
180 Kokkos::parallel_for(
"Utilities_kokkos::GetMatrixDiagonalInverse",
182 KOKKOS_LAMBDA(
const ordinal_type rowIdx) {
183 auto myRow = localMatrix.rowConst(rowIdx);
184 for(ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
185 diagVals(rowIdx, 0) += KAT::magnitude(myRow.value(entryIdx));
187 if(KAT::magnitude(diagVals(rowIdx, 0)) > KAT::magnitude(tol))
188 diagVals(rowIdx, 0) = KAT::one() / diagVals(rowIdx, 0);
190 diagVals(rowIdx, 0) = KAT::zero();
197 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
198 RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
199 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
200 GetMatrixOverlappedDiagonal(
const Matrix& A) {
202 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
203 RCP<Vector> localDiag = VectorFactory::Build(rowMap);
205 const CrsMatrixWrap* crsOp =
dynamic_cast<const CrsMatrixWrap*
>(&A);
207 Teuchos::ArrayRCP<size_t> offsets;
208 crsOp->getLocalDiagOffsets(offsets);
209 crsOp->getLocalDiagCopy(*localDiag,offsets());
212 auto localDiagVals = localDiag->getDeviceLocalView(Xpetra::Access::ReadWrite);
213 const auto diagVals = GetMatrixDiagonal(A)->getDeviceLocalView(Xpetra::Access::ReadOnly);
214 Kokkos::deep_copy(localDiagVals, diagVals);
217 RCP<Vector> diagonal = VectorFactory::Build(colMap);
218 RCP< const Import> importer;
219 importer = A.getCrsGraph()->getImporter();
220 if (importer == Teuchos::null) {
221 importer = ImportFactory::Build(rowMap, colMap);
223 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
228 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
229 void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
230 MyOldScaleMatrix(Matrix& Op,
const Teuchos::ArrayRCP<const SC>& scalingVector,
231 bool doInverse,
bool doFillComplete,
bool doOptimizeStorage)
233 SC one = Teuchos::ScalarTraits<SC>::one();
234 Teuchos::ArrayRCP<SC> sv(scalingVector.size());
236 for (
int i = 0; i < scalingVector.size(); ++i)
237 sv[i] = one / scalingVector[i];
239 for (
int i = 0; i < scalingVector.size(); ++i)
240 sv[i] = scalingVector[i];
243 switch (Op.getRowMap()->lib()) {
244 case Xpetra::UseTpetra:
245 MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
248 case Xpetra::UseEpetra:
251 throw std::runtime_error(
"FIXME");
252 #ifndef __NVCC__ //prevent nvcc warning 257 throw Exceptions::RuntimeError(
"Only Epetra and Tpetra matrices can be scaled.");
258 #ifndef __NVCC__ //prevent nvcc warning 264 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
265 void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& ,
const Teuchos::ArrayRCP<Scalar>& ,
bool ,
bool ) {
266 throw Exceptions::RuntimeError(
"MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
269 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
270 void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Matrix& Op,
const Teuchos::ArrayRCP<SC>& scalingVector,
272 bool doOptimizeStorage)
274 #ifdef HAVE_MUELU_TPETRA 276 Tpetra::CrsMatrix<SC,LO,GO,NO>& tpOp = Op2NonConstTpetraCrs(Op);
278 const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
279 const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
280 const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
282 size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
283 if (maxRowSize == Teuchos::as<size_t>(-1))
286 if (tpOp.isFillComplete())
289 if (Op.isLocallyIndexed() ==
true) {
290 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::local_inds_host_view_type cols;
291 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type vals;
293 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
294 tpOp.getLocalRowView(i, cols, vals);
295 size_t nnz = tpOp.getNumEntriesInLocalRow(i);
296 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::nonconst_values_host_view_type scaledVals(
"scaledVals", nnz);
297 for (
size_t j = 0; j < nnz; ++j)
298 scaledVals[j] = scalingVector[i]*vals[j];
301 tpOp.replaceLocalValues(i, cols, scaledVals);
306 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::global_inds_host_view_type cols;
307 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type vals;
309 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
310 GO gid = rowMap->getGlobalElement(i);
311 tpOp.getGlobalRowView(gid, cols, vals);
312 size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
313 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::nonconst_values_host_view_type scaledVals(
"scaledVals", nnz);
316 for (
size_t j = 0; j < nnz; ++j)
317 scaledVals[j] = scalingVector[i]*vals[j];
320 tpOp.replaceGlobalValues(gid, cols, scaledVals);
325 if (doFillComplete) {
326 if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
327 throw Exceptions::RuntimeError(
"In Utils_kokkos::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
329 RCP<Teuchos::ParameterList> params = rcp(
new Teuchos::ParameterList());
330 params->set(
"Optimize Storage", doOptimizeStorage);
331 params->set(
"No Nonlocal Changes",
true);
332 Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
335 throw Exceptions::RuntimeError(
"Only Tpetra::CrsMatrix types can be scaled (Err.1)");
338 throw Exceptions::RuntimeError(
"Matrix scaling is not possible because Tpetra has not been enabled.");
343 template <
class SC,
class LO,
class GO,
class NO>
346 const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol,
347 const bool count_twos_as_dirichlet) {
348 using impl_scalar_type =
typename Kokkos::ArithTraits<SC>::val_type;
349 using ATS = Kokkos::ArithTraits<impl_scalar_type>;
352 auto localMatrix = A.getLocalMatrixDevice();
353 LO numRows = A.getNodeNumRows();
356 if (count_twos_as_dirichlet)
357 Kokkos::parallel_for(
"MueLu:Utils::DetectDirichletRows_Twos_As_Dirichlet", range_type(0,numRows),
358 KOKKOS_LAMBDA(
const LO row) {
359 auto rowView = localMatrix.row(row);
360 auto length = rowView.length;
362 boundaryNodes(row) =
true;
364 decltype(length) colID;
365 for (colID = 0; colID < length; colID++)
366 if ((rowView.colidx(colID) != row) &&
367 (ATS::magnitude(rowView.value(colID)) > tol)) {
368 if (!boundaryNodes(row))
370 boundaryNodes(row) =
false;
373 boundaryNodes(row) =
true;
377 Kokkos::parallel_for(
"MueLu:Utils::DetectDirichletRows", range_type(0,numRows),
378 KOKKOS_LAMBDA(
const LO row) {
379 auto rowView = localMatrix.row(row);
380 auto length = rowView.length;
382 boundaryNodes(row) =
true;
383 for (decltype(length) colID = 0; colID < length; colID++)
384 if ((rowView.colidx(colID) != row) &&
385 (ATS::magnitude(rowView.value(colID)) > tol)) {
386 boundaryNodes(row) =
false;
391 return boundaryNodes;
394 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
397 DetectDirichletRows(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& tol,
const bool count_twos_as_dirichlet) {
398 return MueLu::DetectDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, tol, count_twos_as_dirichlet);
401 template <
class Node>
404 DetectDirichletRows(
const Xpetra::Matrix<double,int,int,Node>& A,
const typename Teuchos::ScalarTraits<double>::magnitudeType& tol,
const bool count_twos_as_dirichlet) {
405 return MueLu::DetectDirichletRows<double,int,int,Node>(A, tol,count_twos_as_dirichlet);
409 template <
class SC,
class LO,
class GO,
class NO>
413 using ATS = Kokkos::ArithTraits<SC>;
414 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
417 SC zero = ATS::zero();
420 auto localMatrix = A.getLocalMatrixDevice();
421 LO numRows = A.getNodeNumRows();
423 Teuchos::RCP<const Xpetra::Map<LO,GO,NO> > domMap = A.getDomainMap();
424 Teuchos::RCP<const Xpetra::Map<LO,GO,NO> > colMap = A.getColMap();
425 Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > myColsToZero = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(colMap,1);
426 myColsToZero->putScalar(zero);
427 auto myColsToZeroView = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
429 Kokkos::parallel_for(
"MueLu:Utils::DetectDirichletCols1", range_type(0,numRows),
430 KOKKOS_LAMBDA(
const LO row) {
431 if (dirichletRows(row)) {
432 auto rowView = localMatrix.row(row);
433 auto length = rowView.length;
435 for (decltype(length) colID = 0; colID < length; colID++)
436 myColsToZeroView(rowView.colidx(colID),0) = one;
440 Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > globalColsToZero = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(domMap,1);
441 globalColsToZero->putScalar(zero);
442 Teuchos::RCP<Xpetra::Export<LO,GO,NO> > exporter = Xpetra::ExportFactory<LO,GO,NO>::Build(colMap,domMap);
444 globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
446 myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
448 auto myCols = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
449 size_t numColEntries = colMap->getNodeNumElements();
451 const typename ATS::magnitudeType eps = 2.0*ATS::eps();
453 Kokkos::parallel_for(
"MueLu:Utils::DetectDirichletCols2", range_type(0,numColEntries),
454 KOKKOS_LAMBDA (
const size_t i) {
455 dirichletCols(i) = impl_ATS::magnitude(myCols(i,0))>eps;
457 return dirichletCols;
461 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
466 return MueLu::DetectDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
469 template <
class Node>
474 return MueLu::DetectDirichletCols<double,int,int,Node>(A, dirichletRows);
479 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
480 void FindNonZeros(
const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
482 using ATS = Kokkos::ArithTraits<Scalar>;
483 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
485 TEUCHOS_ASSERT(vals.extent(0) == nonzeros.
extent(0));
486 const typename ATS::magnitudeType eps = 2.0*impl_ATS::eps();
488 Kokkos::parallel_for(
"MueLu:Maxwell1::FindNonZeros", range_type(0,vals.extent(0)),
489 KOKKOS_LAMBDA (
const size_t i) {
490 nonzeros(i) = (impl_ATS::magnitude(vals(i,0)) > eps);
494 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
497 FindNonZeros(
const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
499 MueLu::FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(vals,nonzeros);
502 template <
class Node>
505 FindNonZeros(
const typename Xpetra::MultiVector<double,int,int,Node>::dual_view_type::t_dev_const_um vals,
507 MueLu::FindNonZeros<double,int,int,Node>(vals,nonzeros);
511 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
517 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
518 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
519 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
520 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
521 TEUCHOS_ASSERT(dirichletRows.
extent(0) == rowMap->getNodeNumElements());
522 TEUCHOS_ASSERT(dirichletCols.
extent(0) == colMap->getNodeNumElements());
523 TEUCHOS_ASSERT(dirichletDomain.
extent(0) == domMap->getNodeNumElements());
524 RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,
true);
526 auto myColsToZeroView = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
527 auto localMatrix = A.getLocalMatrixDevice();
528 Kokkos::parallel_for(
"MueLu:Maxwell1::DetectDirichletCols", range_type(0,rowMap->getNodeNumElements()),
530 if (dirichletRows(row)) {
531 auto rowView = localMatrix.row(row);
532 auto length = rowView.length;
534 for (decltype(length) colID = 0; colID < length; colID++)
535 myColsToZeroView(rowView.colidx(colID),0) = one;
539 RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
540 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
541 if (!importer.is_null()) {
542 globalColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,
true);
544 globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
546 myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
549 globalColsToZero = myColsToZero;
550 FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(globalColsToZero->getDeviceLocalView(Xpetra::Access::ReadOnly),dirichletDomain);
551 FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(myColsToZero->getDeviceLocalView(Xpetra::Access::ReadOnly),dirichletCols);
555 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
562 MueLu::DetectDirichletColsAndDomains<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A,dirichletRows,dirichletCols,dirichletDomain);
565 template <
class Node>
572 MueLu::DetectDirichletColsAndDomains<double,int,int,Node>(A,dirichletRows,dirichletCols,dirichletDomain);
577 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
584 auto localMatrix = A->getLocalMatrixDevice();
587 Kokkos::parallel_for(
"MueLu:Utils::ZeroDirichletRows", range_type(0,numRows),
589 if (dirichletRows(row)) {
590 auto rowView = localMatrix.row(row);
591 auto length = rowView.length;
592 for (decltype(length) colID = 0; colID < length; colID++)
593 rowView.value(colID) = replaceWith;
598 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
601 ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
604 MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows, replaceWith);
607 template <
class Node>
612 double replaceWith) {
613 return MueLu::ZeroDirichletRows<double,int,int,Node>(A, dirichletRows, replaceWith);
618 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
624 auto myCols = X->getDeviceLocalView(Xpetra::Access::ReadWrite);
625 size_t numVecs = X->getNumVectors();
626 Kokkos::parallel_for(
"MueLu:Utils::ZeroDirichletRows_MV", range_type(0,dirichletRows.size()),
627 KOKKOS_LAMBDA(
const size_t i) {
628 if (dirichletRows(i)) {
629 for(
size_t j=0; j<numVecs; j++)
630 myCols(i,j) = replaceWith;
635 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
638 ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
641 MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X, dirichletRows, replaceWith);
644 template <
class Node>
647 ZeroDirichletRows(RCP<Xpetra::MultiVector<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType,
int,
int,
Node> >& X,
649 double replaceWith) {
650 return MueLu::ZeroDirichletRows<double,int,int,Node>(X, dirichletRows, replaceWith);
655 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
662 auto localMatrix = A->getLocalMatrixDevice();
665 Kokkos::parallel_for(
"MueLu:Utils::ZeroDirichletCols", range_type(0,numRows),
667 auto rowView = localMatrix.row(row);
668 auto length = rowView.length;
669 for (decltype(length) colID = 0; colID < length; colID++)
670 if (dirichletCols(rowView.colidx(colID))) {
671 rowView.value(colID) = replaceWith;
676 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
679 ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
682 MueLu::ZeroDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletCols, replaceWith);
685 template <
class Node>
690 double replaceWith) {
691 return MueLu::ZeroDirichletCols<double,int,int,Node>(A, dirichletCols, replaceWith);
695 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
697 const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
700 typedef Teuchos::ScalarTraits<Scalar> STS;
701 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> rowmap = A.getRowMap();
702 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getNodeNumElements()); ++row) {
703 size_t nnz = A.getNumEntriesInLocalRow(row);
704 ArrayView<const LocalOrdinal> indices;
705 ArrayView<const Scalar> vals;
706 A.getLocalRowView(row, indices, vals);
708 Scalar rowsum = STS::zero();
709 Scalar diagval = STS::zero();
710 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
713 diagval = vals[colID];
714 rowsum += vals[colID];
716 if (STS::real(rowsum) > STS::magnitude(diagval) * rowSumTol)
717 dirichletRows(row) =
true;
721 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
725 const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
728 MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A,rowSumTol,dirichletRows);
732 template <
class Node>
736 const typename Teuchos::ScalarTraits<double>::magnitudeType rowSumTol,
739 MueLu::ApplyRowSumCriterion<double, int, int, Node>(A,rowSumTol,dirichletRows);
743 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
744 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
745 Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
746 RealValuedToScalarMultiVector(RCP<RealValuedMultiVector > X) {
747 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xscalar;
748 #if defined(HAVE_XPETRA_TPETRA) && (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) 752 if ((
typeid(
Scalar).name() ==
typeid(std::complex<double>).name()) ||
753 (
typeid(
Scalar).name() ==
typeid(std::complex<float>).name())) {
754 size_t numVecs = X->getNumVectors();
755 Xscalar = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X->getMap(),numVecs);
756 auto XVec = X->getDeviceLocalView();
757 auto XVecScalar = Xscalar->getDeviceLocalView();
759 Kokkos::parallel_for(
"MueLu:Utils::RealValuedToScalarMultiVector", range_type(0,X->getLocalLength()),
760 KOKKOS_LAMBDA(
const size_t i) {
761 for (
size_t j=0; j<numVecs; j++)
762 XVecScalar(i,j) = XVec(i,j);
766 Xscalar = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(X);
770 template <
class Node>
771 RCP<Xpetra::MultiVector<double,int,int,Node> >
772 Utilities_kokkos<double,int,int,Node>::
773 RealValuedToScalarMultiVector(RCP<Xpetra::MultiVector<Magnitude,int,int,Node> > X) {
778 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
779 Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> >
ReverseCuthillMcKee(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Op) {
780 using local_matrix_type =
typename Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::local_matrix_type;
781 using local_graph_type =
typename local_matrix_type::staticcrsgraph_type;
782 using lno_nnz_view_t =
typename local_graph_type::entries_type::non_const_type;
783 using device =
typename local_graph_type::device_type;
784 using execution_space =
typename local_matrix_type::execution_space;
785 using ordinal_type =
typename local_matrix_type::ordinal_type;
787 local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
789 lno_nnz_view_t rcmOrder = KokkosGraph::Experimental::graph_rcm
790 <device,
typename local_graph_type::row_map_type,
typename local_graph_type::entries_type, lno_nnz_view_t>
791 (localGraph.row_map, localGraph.entries);
793 RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > retval =
794 Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRowMap());
797 auto view1D = Kokkos::subview(retval->getDeviceLocalView(Xpetra::Access::ReadWrite),Kokkos::ALL (), 0);
798 Kokkos::parallel_for(
"Utilities_kokkos::ReverseCuthillMcKee",
800 KOKKOS_LAMBDA(
const ordinal_type rowIdx) {
801 view1D(rcmOrder(rowIdx)) = rowIdx;
806 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
807 Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> >
CuthillMcKee(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Op) {
808 using local_matrix_type =
typename Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::local_matrix_type;
809 using local_graph_type =
typename local_matrix_type::staticcrsgraph_type;
810 using lno_nnz_view_t =
typename local_graph_type::entries_type::non_const_type;
811 using device =
typename local_graph_type::device_type;
812 using execution_space =
typename local_matrix_type::execution_space;
813 using ordinal_type =
typename local_matrix_type::ordinal_type;
815 local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
818 lno_nnz_view_t rcmOrder = KokkosGraph::Experimental::graph_rcm
819 <device,
typename local_graph_type::row_map_type,
typename local_graph_type::entries_type, lno_nnz_view_t>
820 (localGraph.row_map, localGraph.entries);
822 RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > retval =
823 Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRowMap());
826 auto view1D = Kokkos::subview(retval->getDeviceLocalView(Xpetra::Access::ReadWrite),Kokkos::ALL (), 0);
828 Kokkos::parallel_for(
"Utilities_kokkos::ReverseCuthillMcKee",
830 KOKKOS_LAMBDA(
const ordinal_type rowIdx) {
831 view1D(rcmOrder(numRows - 1 - rowIdx)) = rowIdx;
836 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
837 Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> >
839 return MueLu::ReverseCuthillMcKee<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Op);
842 template <
class Node>
843 Teuchos::RCP<Xpetra::Vector<int,int,int,Node> >
845 return MueLu::ReverseCuthillMcKee<double,int,int,Node>(Op);
848 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
849 Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> >
851 return MueLu::CuthillMcKee<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Op);
854 template <
class Node>
855 Teuchos::RCP<Xpetra::Vector<int,int,int,Node> >
857 return MueLu::CuthillMcKee<double,int,int,Node>(Op);
862 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
866 TEUCHOS_ASSERT(A->isFillComplete());
867 using ATS = Kokkos::ArithTraits<Scalar>;
868 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
871 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A->getDomainMap();
872 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ranMap = A->getRangeMap();
873 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Rmap = A->getRowMap();
874 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Cmap = A->getColMap();
876 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getNodeNumElements());
878 const Scalar one = impl_ATS::one();
879 const Scalar zero = impl_ATS::zero();
881 auto localMatrix = A->getLocalMatrixDevice();
882 auto localRmap = Rmap->getLocalMap();
883 auto localCmap = Cmap->getLocalMap();
885 Kokkos::parallel_for(
"MueLu::Utils::ApplyOAZ",range_type(0,dirichletRows.
extent(0)),
887 if (dirichletRows(row)){
888 auto rowView = localMatrix.row(row);
889 auto length = rowView.length;
890 auto row_gid = localRmap.getGlobalElement(row);
891 auto row_lid = localCmap.getLocalElement(row_gid);
893 for (decltype(length) colID = 0; colID < length; colID++)
894 if (rowView.colidx(colID) == row_lid)
895 rowView.value(colID) = one;
897 rowView.value(colID) = zero;
902 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
907 MueLu::ApplyOAZToMatrixRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
910 template <
class Node>
915 MueLu::ApplyOAZToMatrixRows<double,int,int,Node>(A, dirichletRows);
920 #define MUELU_UTILITIES_KOKKOS_SHORT 921 #endif // MUELU_UTILITIES_KOKKOS_DEF_HPP MueLu::DefaultLocalOrdinal LocalOrdinal
Teuchos::RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< std::is_integral< iType >::value, size_t >::type extent(const iType &r) const noexcept
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
void ZeroDirichletRows(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, Scalar replaceWith)
Teuchos::RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Namespace for MueLu classes and methods.
void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows)
void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Kokkos::View< bool *, typename Node::device_type > &dirichletRows, Kokkos::View< bool *, typename Node::device_type > dirichletCols, Kokkos::View< bool *, typename Node::device_type > dirichletDomain)
void ZeroDirichletRows(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &X, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, Scalar replaceWith)
Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
MueLu::DefaultScalar Scalar
void FindNonZeros(const typename Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >::dual_view_type::t_dev_const_um vals, Kokkos::View< bool *, typename Node::device_type > nonzeros)
void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, typename Node::device_type > &dirichletRows)
Kokkos::View< bool *, typename NO::device_type > DetectDirichletCols(const Xpetra::Matrix< SC, LO, GO, NO > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows)