MueLu  Version of the Day
MueLu_Utilities_kokkos_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_UTILITIES_KOKKOS_DEF_HPP
47 #define MUELU_UTILITIES_KOKKOS_DEF_HPP
48 
49 #include <algorithm>
50 #include <Teuchos_DefaultComm.hpp>
51 #include <Teuchos_ParameterList.hpp>
52 
53 #include "MueLu_ConfigDefs.hpp"
54 
55 #ifdef HAVE_MUELU_EPETRA
56 # ifdef HAVE_MPI
57 # include "Epetra_MpiComm.h"
58 # endif
59 #endif
60 
61 #include <Kokkos_Core.hpp>
62 #include <KokkosSparse_CrsMatrix.hpp>
63 #include <KokkosSparse_getDiagCopy.hpp>
64 
65 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
66 #include <EpetraExt_MatrixMatrix.h>
67 #include <EpetraExt_RowMatrixOut.h>
69 #include <EpetraExt_CrsMatrixIn.h>
71 #include <EpetraExt_BlockMapIn.h>
72 #include <Xpetra_EpetraUtils.hpp>
73 #include <Xpetra_EpetraMultiVector.hpp>
74 #include <EpetraExt_BlockMapOut.h>
75 #endif
76 
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>
84 #endif
85 
86 #ifdef HAVE_MUELU_EPETRA
87 #include <Xpetra_EpetraMap.hpp>
88 #endif
89 
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>
104 
106 
107 #include <KokkosKernels_Handle.hpp>
108 #include <KokkosGraph_RCM.hpp>
109 
110 
111 namespace MueLu {
112 
113 
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);
120 
121  A.getLocalDiagCopy(*diag);
122 
123  return diag;
124  } //GetMatrixDiagonal
125 
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");
131  // Some useful type definitions
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;
138  // Be careful with this one, if using Kokkos::ArithTraits<Scalar>
139  // you are likely to run into errors when handling std::complex<>
140  // a good way to work around that is to use the following:
141  // using KAT = Kokkos::ArithTraits<Kokkos::ArithTraits<Scalar>::val_type> >
142  // here we have: value_type = Kokkos::ArithTraits<Scalar>::val_type
143  using KAT = Kokkos::ArithTraits<value_type>;
144 
145  // Get/Create distributed objects
146  RCP<const Map> rowMap = A.getRowMap();
147  RCP<Vector> diag = VectorFactory::Build(rowMap,false);
148 
149  // Now generate local objects
150  local_matrix_type localMatrix = A.getLocalMatrixDevice();
151  auto diagVals = diag->getDeviceLocalView(Xpetra::Access::ReadWrite);
152 
153  ordinal_type numRows = localMatrix.graph.numRows();
154 
155  // Note: 2019-11-21, LBV
156  // This could be implemented with a TeamPolicy over the rows
157  // and a TeamVectorRange over the entries in a row if performance
158  // becomes more important here.
159  if (!doLumped)
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);
170  } else {
171  diagVals(rowIdx, 0) = KAT::zero();
172  }
173  break;
174  }
175  }
176 
177  if(!foundDiagEntry) {diagVals(rowIdx, 0) = KAT::zero();}
178  });
179  else
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));
186  }
187  if(KAT::magnitude(diagVals(rowIdx, 0)) > KAT::magnitude(tol))
188  diagVals(rowIdx, 0) = KAT::one() / diagVals(rowIdx, 0);
189  else
190  diagVals(rowIdx, 0) = KAT::zero();
191 
192  });
193 
194  return diag;
195  } //GetMatrixDiagonalInverse
196 
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) {
201  // FIXME_KOKKOS
202  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
203  RCP<Vector> localDiag = VectorFactory::Build(rowMap);
204 
205  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
206  if (crsOp != NULL) {
207  Teuchos::ArrayRCP<size_t> offsets;
208  crsOp->getLocalDiagOffsets(offsets);
209  crsOp->getLocalDiagCopy(*localDiag,offsets());
210  }
211  else {
212  auto localDiagVals = localDiag->getDeviceLocalView(Xpetra::Access::ReadWrite);
213  const auto diagVals = GetMatrixDiagonal(A)->getDeviceLocalView(Xpetra::Access::ReadOnly);
214  Kokkos::deep_copy(localDiagVals, diagVals);
215  }
216 
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);
222  }
223  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
224 
225  return diagonal;
226  } //GetMatrixOverlappedDiagonal
227 
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)
232  {
233  SC one = Teuchos::ScalarTraits<SC>::one();
234  Teuchos::ArrayRCP<SC> sv(scalingVector.size());
235  if (doInverse) {
236  for (int i = 0; i < scalingVector.size(); ++i)
237  sv[i] = one / scalingVector[i];
238  } else {
239  for (int i = 0; i < scalingVector.size(); ++i)
240  sv[i] = scalingVector[i];
241  }
242 
243  switch (Op.getRowMap()->lib()) {
244  case Xpetra::UseTpetra:
245  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
246  break;
247 
248  case Xpetra::UseEpetra:
249  // FIXME_KOKKOS
250  // Utils2_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
251  throw std::runtime_error("FIXME");
252 #ifndef __NVCC__ //prevent nvcc warning
253  break;
254 #endif
255 
256  default:
257  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
258 #ifndef __NVCC__ //prevent nvcc warning
259  break;
260 #endif
261  }
262  }
263 
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>& /* Op */, const Teuchos::ArrayRCP<Scalar>& /* scalingVector */, bool /* doFillComplete */, bool /* doOptimizeStorage */) {
266  throw Exceptions::RuntimeError("MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
267  }
268 
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,
271  bool doFillComplete,
272  bool doOptimizeStorage)
273  {
274 #ifdef HAVE_MUELU_TPETRA
275  try {
276  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpOp = Op2NonConstTpetraCrs(Op);
277 
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();
281 
282  size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
283  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
284  maxRowSize = 20;
285 
286  if (tpOp.isFillComplete())
287  tpOp.resumeFill();
288 
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;
292 
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];
299 
300  if (nnz > 0) {
301  tpOp.replaceLocalValues(i, cols, scaledVals);
302  }
303  } //for (size_t i=0; ...
304 
305  } else {
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;
308 
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);
314 
315  // FIXME FIXME FIXME FIXME FIXME FIXME
316  for (size_t j = 0; j < nnz; ++j)
317  scaledVals[j] = scalingVector[i]*vals[j]; //FIXME i or gid?
318 
319  if (nnz > 0) {
320  tpOp.replaceGlobalValues(gid, cols, scaledVals);
321  }
322  } //for (size_t i=0; ...
323  }
324 
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");
328 
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);
333  }
334  } catch(...) {
335  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
336  }
337 #else
338  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
339 #endif
340  } //MyOldScaleMatrix_Tpetra()
341 
342 
343  template <class SC, class LO, class GO, class NO>
345  DetectDirichletRows(const Xpetra::Matrix<SC,LO,GO,NO>& A,
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>;
351 
352  auto localMatrix = A.getLocalMatrixDevice();
353  LO numRows = A.getNodeNumRows();
354 
355  Kokkos::View<bool*, typename NO::device_type> boundaryNodes(Kokkos::ViewAllocateWithoutInitializing("boundaryNodes"), numRows);
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;
361 
362  boundaryNodes(row) = true;
363  if (length > 2) {
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))
369  break;
370  boundaryNodes(row) = false;
371  }
372  if (colID == length)
373  boundaryNodes(row) = true;
374  }
375  });
376  else
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;
381 
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;
387  break;
388  }
389  });
390 
391  return boundaryNodes;
392  }
393 
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);
399  }
400 
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);
406  }
407 
408 
409  template <class SC, class LO, class GO, class NO>
411  DetectDirichletCols(const Xpetra::Matrix<SC,LO,GO,NO>& A,
413  using ATS = Kokkos::ArithTraits<SC>;
414  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
416 
417  SC zero = ATS::zero();
418  SC one = ATS::one();
419 
420  auto localMatrix = A.getLocalMatrixDevice();
421  LO numRows = A.getNodeNumRows();
422 
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);
428  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
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;
434 
435  for (decltype(length) colID = 0; colID < length; colID++)
436  myColsToZeroView(rowView.colidx(colID),0) = one;
437  }
438  });
439 
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);
443  // export to domain map
444  globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
445  // import to column map
446  myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
447 
448  auto myCols = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
449  size_t numColEntries = colMap->getNodeNumElements();
450  Kokkos::View<bool*, typename NO::device_type> dirichletCols(Kokkos::ViewAllocateWithoutInitializing("dirichletCols"), numColEntries);
451  const typename ATS::magnitudeType eps = 2.0*ATS::eps();
452 
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;
456  });
457  return dirichletCols;
458  }
459 
460 
461  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
464  DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
466  return MueLu::DetectDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
467  }
468 
469  template <class Node>
472  DetectDirichletCols(const Xpetra::Matrix<double,int,int,Node>& A,
474  return MueLu::DetectDirichletCols<double,int,int,Node>(A, dirichletRows);
475  }
476 
477 
478  // Find Nonzeros in a device view
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();
487 
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);
491  });
492  }
493 
494  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
495  void
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);
500  }
501 
502  template <class Node>
503  void
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);
508  }
509 
510  // Detect Dirichlet cols and domains
511  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
512  void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
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, /*zeroOut=*/true);
525  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
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()),
529  KOKKOS_LAMBDA(const LocalOrdinal row) {
530  if (dirichletRows(row)) {
531  auto rowView = localMatrix.row(row);
532  auto length = rowView.length;
533 
534  for (decltype(length) colID = 0; colID < length; colID++)
535  myColsToZeroView(rowView.colidx(colID),0) = one;
536  }
537  });
538 
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, /*zeroOut=*/true);
543  // export to domain map
544  globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
545  // import to column map
546  myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
547  }
548  else
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);
552  }
553 
554 
555  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
556  void
558  DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
562  MueLu::DetectDirichletColsAndDomains<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A,dirichletRows,dirichletCols,dirichletDomain);
563  }
564 
565  template <class Node>
566  void
568  DetectDirichletColsAndDomains(const Xpetra::Matrix<double,int,int,Node>& A,
572  MueLu::DetectDirichletColsAndDomains<double,int,int,Node>(A,dirichletRows,dirichletCols,dirichletDomain);
573  }
574 
575 
576  // Zeros out rows
577  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
578  void
579  ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& A,
581  Scalar replaceWith) {
583 
584  auto localMatrix = A->getLocalMatrixDevice();
585  LocalOrdinal numRows = A->getNodeNumRows();
586 
587  Kokkos::parallel_for("MueLu:Utils::ZeroDirichletRows", range_type(0,numRows),
588  KOKKOS_LAMBDA(const LocalOrdinal row) {
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;
594  }
595  });
596  }
597 
598  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
599  void
601  ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
603  Scalar replaceWith) {
604  MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows, replaceWith);
605  }
606 
607  template <class Node>
608  void
610  ZeroDirichletRows(RCP<Xpetra::Matrix<double, int, int, Node> >& A,
612  double replaceWith) {
613  return MueLu::ZeroDirichletRows<double,int,int,Node>(A, dirichletRows, replaceWith);
614  }
615 
616 
617  // Zeros out rows
618  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
619  void
620  ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
622  Scalar replaceWith) {
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;
631  }
632  });
633  }
634 
635  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
636  void
638  ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
640  Scalar replaceWith) {
641  MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X, dirichletRows, replaceWith);
642  }
643 
644  template <class Node>
645  void
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);
651  }
652 
653 
654  // Zeros out columns
655  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
656  void
657  ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
659  Scalar replaceWith) {
661 
662  auto localMatrix = A->getLocalMatrixDevice();
663  LocalOrdinal numRows = A->getNodeNumRows();
664 
665  Kokkos::parallel_for("MueLu:Utils::ZeroDirichletCols", range_type(0,numRows),
666  KOKKOS_LAMBDA(const LocalOrdinal row) {
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;
672  }
673  });
674  }
675 
676  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
677  void
679  ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
681  Scalar replaceWith) {
682  MueLu::ZeroDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletCols, replaceWith);
683  }
684 
685  template <class Node>
686  void
688  ZeroDirichletCols(RCP<Xpetra::Matrix<double,int,int,Node> >& A,
690  double replaceWith) {
691  return MueLu::ZeroDirichletCols<double,int,int,Node>(A, dirichletCols, replaceWith);
692  }
693 
694  // Applies rowsum criterion
695  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
696  void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
697  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
699  {
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);
707 
708  Scalar rowsum = STS::zero();
709  Scalar diagval = STS::zero();
710  for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
711  LocalOrdinal col = indices[colID];
712  if (row == col)
713  diagval = vals[colID];
714  rowsum += vals[colID];
715  }
716  if (STS::real(rowsum) > STS::magnitude(diagval) * rowSumTol)
717  dirichletRows(row) = true;
718  }
719  }
720 
721  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
722  void
724  ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
725  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
727  {
728  MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A,rowSumTol,dirichletRows);
729  }
730 
731 
732  template <class Node>
733  void
735  ApplyRowSumCriterion(const Xpetra::Matrix<double,int,int,Node>& A,
736  const typename Teuchos::ScalarTraits<double>::magnitudeType rowSumTol,
738  {
739  MueLu::ApplyRowSumCriterion<double, int, int, Node>(A,rowSumTol,dirichletRows);
740  }
741 
742 
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))
750 
751  // Need to cast the real-valued multivector to Scalar=complex
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();
758 
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);
763  });
764  } else
765 #endif
766  Xscalar = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(X);
767  return Xscalar;
768  }
769 
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) {
774  return X;
775  }
776 
777 
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;
786 
787  local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
788 
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);
792 
793  RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > retval =
794  Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRowMap());
795 
796  // Copy out and reorder data
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;
802  });
803  return retval;
804  }
805 
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;
814 
815  local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
816  LocalOrdinal numRows = localGraph.numRows();
817 
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);
821 
822  RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > retval =
823  Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRowMap());
824 
825  // Copy out data
826  auto view1D = Kokkos::subview(retval->getDeviceLocalView(Xpetra::Access::ReadWrite),Kokkos::ALL (), 0);
827  // Since KokkosKernels produced RCM, also reverse the order of the view to get CM
828  Kokkos::parallel_for("Utilities_kokkos::ReverseCuthillMcKee",
830  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
831  view1D(rcmOrder(numRows - 1 - rowIdx)) = rowIdx;
832  });
833  return retval;
834  }
835 
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);
840  }
841 
842  template <class Node>
843  Teuchos::RCP<Xpetra::Vector<int,int,int,Node> >
845  return MueLu::ReverseCuthillMcKee<double,int,int,Node>(Op);
846  }
847 
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);
852  }
853 
854  template <class Node>
855  Teuchos::RCP<Xpetra::Vector<int,int,int,Node> >
857  return MueLu::CuthillMcKee<double,int,int,Node>(Op);
858  }
859 
860  // Applies Ones-and-Zeros to matrix rows
861  // Takes a Boolean array.
862  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
863  void
864  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
866  TEUCHOS_ASSERT(A->isFillComplete());
867  using ATS = Kokkos::ArithTraits<Scalar>;
868  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
870 
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();
875 
876  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getNodeNumElements());
877 
878  const Scalar one = impl_ATS::one();
879  const Scalar zero = impl_ATS::zero();
880 
881  auto localMatrix = A->getLocalMatrixDevice();
882  auto localRmap = Rmap->getLocalMap();
883  auto localCmap = Cmap->getLocalMap();
884 
885  Kokkos::parallel_for("MueLu::Utils::ApplyOAZ",range_type(0,dirichletRows.extent(0)),
886  KOKKOS_LAMBDA(const LocalOrdinal row) {
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);
892 
893  for (decltype(length) colID = 0; colID < length; colID++)
894  if (rowView.colidx(colID) == row_lid)
895  rowView.value(colID) = one;
896  else
897  rowView.value(colID) = zero;
898  }
899  });
900  }
901 
902  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
903  void
905  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
907  MueLu::ApplyOAZToMatrixRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
908  }
909 
910  template <class Node>
911  void
913  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<double,int,int,Node> >& A,
915  MueLu::ApplyOAZToMatrixRows<double,int,int,Node>(A, dirichletRows);
916  }
917 
918 } //namespace MueLu
919 
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)
MueLu::DefaultNode Node
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)