MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
68 #include <EpetraExt_MultiVectorOut.h>
69 #include <EpetraExt_CrsMatrixIn.h>
70 #include <EpetraExt_MultiVectorIn.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 
105 #include <MueLu_Utilities_kokkos_decl.hpp>
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  GetMatrixDiagonalInverse(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
129  typename Teuchos::ScalarTraits<Scalar>::magnitudeType tol, const bool doLumped) {
130  Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer("Utilities_kokkos::GetMatrixDiagonalInverse");
131  // Some useful type definitions
132  using Matrix = Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
133  using Map = Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>;
134  using Vector = Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
135  using VectorFactory = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
136  using local_matrix_type = typename Matrix::local_matrix_type;
137  // using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
138  using value_type = typename local_matrix_type::value_type;
139  using ordinal_type = typename local_matrix_type::ordinal_type;
140  using execution_space = typename local_matrix_type::execution_space;
141  // using memory_space = typename local_matrix_type::memory_space;
142  // Be careful with this one, if using Kokkos::ArithTraits<Scalar>
143  // you are likely to run into errors when handling std::complex<>
144  // a good way to work around that is to use the following:
145  // using KAT = Kokkos::ArithTraits<Kokkos::ArithTraits<Scalar>::val_type> >
146  // here we have: value_type = Kokkos::ArithTraits<Scalar>::val_type
147  using KAT = Kokkos::ArithTraits<value_type>;
148 
149  // Get/Create distributed objects
150  RCP<const Map> rowMap = A.getRowMap();
151  RCP<Vector> diag = VectorFactory::Build(rowMap,false);
152 
153  // Now generate local objects
154  local_matrix_type localMatrix = A.getLocalMatrixDevice();
155  auto diagVals = diag->getDeviceLocalView(Xpetra::Access::ReadWrite);
156 
157  ordinal_type numRows = localMatrix.graph.numRows();
158 
159  // Note: 2019-11-21, LBV
160  // This could be implemented with a TeamPolicy over the rows
161  // and a TeamVectorRange over the entries in a row if performance
162  // becomes more important here.
163  if (!doLumped)
164  Kokkos::parallel_for("Utilities_kokkos::GetMatrixDiagonalInverse",
165  Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
166  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
167  bool foundDiagEntry = false;
168  auto myRow = localMatrix.rowConst(rowIdx);
169  for(ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
170  if(myRow.colidx(entryIdx) == rowIdx) {
171  foundDiagEntry = true;
172  if(KAT::magnitude(myRow.value(entryIdx)) > KAT::magnitude(tol)) {
173  diagVals(rowIdx, 0) = KAT::one() / myRow.value(entryIdx);
174  } else {
175  diagVals(rowIdx, 0) = KAT::zero();
176  }
177  break;
178  }
179  }
180 
181  if(!foundDiagEntry) {diagVals(rowIdx, 0) = KAT::zero();}
182  });
183  else
184  Kokkos::parallel_for("Utilities_kokkos::GetMatrixDiagonalInverse",
185  Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
186  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
187  auto myRow = localMatrix.rowConst(rowIdx);
188  for(ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
189  diagVals(rowIdx, 0) += KAT::magnitude(myRow.value(entryIdx));
190  }
191  if(KAT::magnitude(diagVals(rowIdx, 0)) > KAT::magnitude(tol))
192  diagVals(rowIdx, 0) = KAT::one() / diagVals(rowIdx, 0);
193  else
194  diagVals(rowIdx, 0) = KAT::zero();
195 
196  });
197 
198  return diag;
199  } //GetMatrixDiagonalInverse
200 
201  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
202  Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal,Node> >
203  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal,Node>::
204  GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol, const bool doLumped)
205  {
206  return MueLu::GetMatrixDiagonalInverse<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A,tol,doLumped);
207  }
208 
209  template <class Node>
210  Teuchos::RCP<Xpetra::Vector<double,int,int,Node> >
211  Utilities_kokkos<double,int,int,Node>::
212  GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol, const bool doLumped)
213  {
214  return MueLu::GetMatrixDiagonalInverse<double, int, int, Node>(A,tol,doLumped);
215  }
216 
217  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
218  RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
219  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
220  GetMatrixOverlappedDiagonal(const Matrix& A) {
221  // FIXME_KOKKOS
222  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
223  RCP<Vector> localDiag = VectorFactory::Build(rowMap);
224 
225  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
226  if (crsOp != NULL) {
227  Teuchos::ArrayRCP<size_t> offsets;
228  crsOp->getLocalDiagOffsets(offsets);
229  crsOp->getLocalDiagCopy(*localDiag,offsets());
230  }
231  else {
232  auto localDiagVals = localDiag->getDeviceLocalView(Xpetra::Access::ReadWrite);
233  const auto diagVals = GetMatrixDiagonal(A)->getDeviceLocalView(Xpetra::Access::ReadOnly);
234  Kokkos::deep_copy(localDiagVals, diagVals);
235  }
236 
237  RCP<Vector> diagonal = VectorFactory::Build(colMap);
238  RCP< const Import> importer;
239  importer = A.getCrsGraph()->getImporter();
240  if (importer == Teuchos::null) {
241  importer = ImportFactory::Build(rowMap, colMap);
242  }
243  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
244 
245  return diagonal;
246  } //GetMatrixOverlappedDiagonal
247 
248  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
249  void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
250  MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector,
251  bool doInverse, bool doFillComplete, bool doOptimizeStorage)
252  {
253  SC one = Teuchos::ScalarTraits<SC>::one();
254  Teuchos::ArrayRCP<SC> sv(scalingVector.size());
255  if (doInverse) {
256  for (int i = 0; i < scalingVector.size(); ++i)
257  sv[i] = one / scalingVector[i];
258  } else {
259  for (int i = 0; i < scalingVector.size(); ++i)
260  sv[i] = scalingVector[i];
261  }
262 
263  switch (Op.getRowMap()->lib()) {
264  case Xpetra::UseTpetra:
265  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
266  break;
267 
268  case Xpetra::UseEpetra:
269  // FIXME_KOKKOS
270  // Utils2_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
271  throw std::runtime_error("FIXME");
272 #ifndef __NVCC__ //prevent nvcc warning
273  break;
274 #endif
275 
276  default:
277  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
278 #ifndef __NVCC__ //prevent nvcc warning
279  break;
280 #endif
281  }
282  }
283 
284  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
285  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 */) {
286  throw Exceptions::RuntimeError("MyOldScaleMatrix_Epetra: Epetra needs SC=double and LO=GO=int.");
287  }
288 
289  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
290  void Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
291  bool doFillComplete,
292  bool doOptimizeStorage)
293  {
294 #ifdef HAVE_MUELU_TPETRA
295  try {
296  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpOp = Op2NonConstTpetraCrs(Op);
297 
298  const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
299  const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
300  const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
301 
302  size_t maxRowSize = tpOp.getLocalMaxNumRowEntries();
303  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
304  maxRowSize = 20;
305 
306  if (tpOp.isFillComplete())
307  tpOp.resumeFill();
308 
309  if (Op.isLocallyIndexed() == true) {
310  typename Tpetra::CrsMatrix<SC,LO,GO,NO>::local_inds_host_view_type cols;
311  typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type vals;
312 
313  for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
314  tpOp.getLocalRowView(i, cols, vals);
315  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
316  typename Tpetra::CrsMatrix<SC,LO,GO,NO>::nonconst_values_host_view_type scaledVals("scaledVals", nnz);
317  for (size_t j = 0; j < nnz; ++j)
318  scaledVals[j] = scalingVector[i]*vals[j];
319 
320  if (nnz > 0) {
321  tpOp.replaceLocalValues(i, cols, scaledVals);
322  }
323  } //for (size_t i=0; ...
324 
325  } else {
326  typename Tpetra::CrsMatrix<SC,LO,GO,NO>::global_inds_host_view_type cols;
327  typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type vals;
328 
329  for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
330  GO gid = rowMap->getGlobalElement(i);
331  tpOp.getGlobalRowView(gid, cols, vals);
332  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
333  typename Tpetra::CrsMatrix<SC,LO,GO,NO>::nonconst_values_host_view_type scaledVals("scaledVals", nnz);
334 
335  // FIXME FIXME FIXME FIXME FIXME FIXME
336  for (size_t j = 0; j < nnz; ++j)
337  scaledVals[j] = scalingVector[i]*vals[j]; //FIXME i or gid?
338 
339  if (nnz > 0) {
340  tpOp.replaceGlobalValues(gid, cols, scaledVals);
341  }
342  } //for (size_t i=0; ...
343  }
344 
345  if (doFillComplete) {
346  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
347  throw Exceptions::RuntimeError("In Utils_kokkos::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
348 
349  RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
350  params->set("Optimize Storage", doOptimizeStorage);
351  params->set("No Nonlocal Changes", true);
352  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
353  }
354  } catch(...) {
355  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
356  }
357 #else
358  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
359 #endif
360  } //MyOldScaleMatrix_Tpetra()
361 
362 
363  template <class SC, class LO, class GO, class NO>
364  Kokkos::View<bool*, typename NO::device_type>
365  DetectDirichletRows(const Xpetra::Matrix<SC,LO,GO,NO>& A,
366  const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol,
367  const bool count_twos_as_dirichlet) {
368  using impl_scalar_type = typename Kokkos::ArithTraits<SC>::val_type;
369  using ATS = Kokkos::ArithTraits<impl_scalar_type>;
370  using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
371  using helpers = Xpetra::Helpers<SC,LO,GO,NO>;
372 
373 
374  if(helpers::isTpetraBlockCrs(A)) {
375 #ifdef HAVE_MUELU_TPETRA
376  const Tpetra::BlockCrsMatrix<SC,LO,GO,NO> & Am = helpers::Op2TpetraBlockCrs(A);
377  auto b_graph = Am.getCrsGraph().getLocalGraphDevice();
378  auto b_rowptr = Am.getCrsGraph().getLocalRowPtrsDevice();
379  auto values = Am.getValuesDevice();
380  LO numBlockRows = Am.getLocalNumRows();
381  const LO stride = Am.getBlockSize() * Am.getBlockSize();
382 
383  Kokkos::View<bool*, typename NO::device_type> boundaryNodes(Kokkos::ViewAllocateWithoutInitializing("boundaryNodes"), numBlockRows);
384 
385  if (count_twos_as_dirichlet)
386  throw Exceptions::RuntimeError("BlockCrs does not support counting twos as Dirichlet");
387 
388  Kokkos::parallel_for("MueLu:Utils::DetectDirichletRowsBlockCrs", range_type(0,numBlockRows),
389  KOKKOS_LAMBDA(const LO row) {
390  auto rowView = b_graph.rowConst(row);
391  auto length = rowView.length;
392  LO valstart = b_rowptr[row] * stride;
393 
394  boundaryNodes(row) = true;
395  decltype(length) colID =0;
396  for (; colID < length; colID++) {
397  if (rowView.colidx(colID) != row) {
398  LO current = valstart + colID*stride;
399  for(LO k=0; k<stride; k++) {
400  if (ATS::magnitude(values[current+ k]) > tol) {
401  boundaryNodes(row) = false;
402  break;
403  }
404  }
405  }
406  if(boundaryNodes(row) == false)
407  break;
408  }
409  });
410 
411  return boundaryNodes;
412 #else
413  throw Exceptions::RuntimeError("BlockCrs requires Tpetra");
414 #endif
415  }
416  else {
417  auto localMatrix = A.getLocalMatrixDevice();
418  LO numRows = A.getLocalNumRows();
419  Kokkos::View<bool*, typename NO::device_type> boundaryNodes(Kokkos::ViewAllocateWithoutInitializing("boundaryNodes"), numRows);
420 
421  if (count_twos_as_dirichlet)
422  Kokkos::parallel_for("MueLu:Utils::DetectDirichletRows_Twos_As_Dirichlet", range_type(0,numRows),
423  KOKKOS_LAMBDA(const LO row) {
424  auto rowView = localMatrix.row(row);
425  auto length = rowView.length;
426 
427  boundaryNodes(row) = true;
428  if (length > 2) {
429  decltype(length) colID =0;
430  for ( ; colID < length; colID++)
431  if ((rowView.colidx(colID) != row) &&
432  (ATS::magnitude(rowView.value(colID)) > tol)) {
433  if (!boundaryNodes(row))
434  break;
435  boundaryNodes(row) = false;
436  }
437  if (colID == length)
438  boundaryNodes(row) = true;
439  }
440  });
441  else
442  Kokkos::parallel_for("MueLu:Utils::DetectDirichletRows", range_type(0,numRows),
443  KOKKOS_LAMBDA(const LO row) {
444  auto rowView = localMatrix.row(row);
445  auto length = rowView.length;
446 
447  boundaryNodes(row) = true;
448  for (decltype(length) colID = 0; colID < length; colID++)
449  if ((rowView.colidx(colID) != row) &&
450  (ATS::magnitude(rowView.value(colID)) > tol)) {
451  boundaryNodes(row) = false;
452  break;
453  }
454  });
455  return boundaryNodes;
456  }
457 
458  }
459 
460  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
461  Kokkos::View<bool*, typename Node::device_type>
462  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
463  DetectDirichletRows(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& tol, const bool count_twos_as_dirichlet) {
464  return MueLu::DetectDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, tol, count_twos_as_dirichlet);
465  }
466 
467  template <class Node>
468  Kokkos::View<bool*, typename Node::device_type>
469  Utilities_kokkos<double,int,int,Node>::
470  DetectDirichletRows(const Xpetra::Matrix<double,int,int,Node>& A, const typename Teuchos::ScalarTraits<double>::magnitudeType& tol, const bool count_twos_as_dirichlet) {
471  return MueLu::DetectDirichletRows<double,int,int,Node>(A, tol,count_twos_as_dirichlet);
472  }
473 
474 
475  template <class SC, class LO, class GO, class NO>
476  Kokkos::View<bool*, typename NO::device_type>
477  DetectDirichletCols(const Xpetra::Matrix<SC,LO,GO,NO>& A,
478  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows) {
479  using ATS = Kokkos::ArithTraits<SC>;
480  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
481  using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
482 
483  SC zero = ATS::zero();
484 
485  auto localMatrix = A.getLocalMatrixDevice();
486  LO numRows = A.getLocalNumRows();
487 
488  Teuchos::RCP<const Xpetra::Map<LO,GO,NO> > domMap = A.getDomainMap();
489  Teuchos::RCP<const Xpetra::Map<LO,GO,NO> > colMap = A.getColMap();
490  Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > myColsToZero = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(colMap,1);
491  myColsToZero->putScalar(zero);
492  auto myColsToZeroView = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
493  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
494  Kokkos::parallel_for("MueLu:Utils::DetectDirichletCols1", range_type(0,numRows),
495  KOKKOS_LAMBDA(const LO row) {
496  if (dirichletRows(row)) {
497  auto rowView = localMatrix.row(row);
498  auto length = rowView.length;
499 
500  for (decltype(length) colID = 0; colID < length; colID++)
501  myColsToZeroView(rowView.colidx(colID),0) = impl_ATS::one();
502  }
503  });
504 
505  Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > globalColsToZero = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(domMap,1);
506  globalColsToZero->putScalar(zero);
507  Teuchos::RCP<Xpetra::Export<LO,GO,NO> > exporter = Xpetra::ExportFactory<LO,GO,NO>::Build(colMap,domMap);
508  // export to domain map
509  globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
510  // import to column map
511  myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
512 
513  auto myCols = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadOnly);
514  size_t numColEntries = colMap->getLocalNumElements();
515  Kokkos::View<bool*, typename NO::device_type> dirichletCols(Kokkos::ViewAllocateWithoutInitializing("dirichletCols"), numColEntries);
516  const typename ATS::magnitudeType eps = 2.0*ATS::eps();
517 
518  Kokkos::parallel_for("MueLu:Utils::DetectDirichletCols2", range_type(0,numColEntries),
519  KOKKOS_LAMBDA (const size_t i) {
520  dirichletCols(i) = impl_ATS::magnitude(myCols(i,0))>eps;
521  });
522  return dirichletCols;
523  }
524 
525 
526  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
527  Kokkos::View<bool*, typename Node::device_type>
528  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
529  DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
530  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows) {
531  return MueLu::DetectDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
532  }
533 
534  template <class Node>
535  Kokkos::View<bool*, typename Node::device_type>
536  Utilities_kokkos<double,int,int,Node>::
537  DetectDirichletCols(const Xpetra::Matrix<double,int,int,Node>& A,
538  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
539  return MueLu::DetectDirichletCols<double,int,int,Node>(A, dirichletRows);
540  }
541 
542 
543  // Find Nonzeros in a device view
544  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
545  void FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
546  Kokkos::View<bool*, typename Node::device_type> nonzeros) {
547  using ATS = Kokkos::ArithTraits<Scalar>;
548  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
549  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
550  TEUCHOS_ASSERT(vals.extent(0) == nonzeros.extent(0));
551  const typename ATS::magnitudeType eps = 2.0*impl_ATS::eps();
552 
553  Kokkos::parallel_for("MueLu:Maxwell1::FindNonZeros", range_type(0,vals.extent(0)),
554  KOKKOS_LAMBDA (const size_t i) {
555  nonzeros(i) = (impl_ATS::magnitude(vals(i,0)) > eps);
556  });
557  }
558 
559  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
560  void
561  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
562  FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
563  Kokkos::View<bool*, typename Node::device_type> nonzeros) {
564  MueLu::FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(vals,nonzeros);
565  }
566 
567  template <class Node>
568  void
569  Utilities_kokkos<double,int,int,Node>::
570  FindNonZeros(const typename Xpetra::MultiVector<double,int,int,Node>::dual_view_type::t_dev_const_um vals,
571  Kokkos::View<bool*, typename Node::device_type> nonzeros) {
572  MueLu::FindNonZeros<double,int,int,Node>(vals,nonzeros);
573  }
574 
575  // Detect Dirichlet cols and domains
576  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
577  void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
578  const Kokkos::View<bool*, typename Node::device_type> & dirichletRows,
579  Kokkos::View<bool*, typename Node::device_type> dirichletCols,
580  Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
581  using ATS = Kokkos::ArithTraits<Scalar>;
582  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
583  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
584  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
585  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
586  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
587  TEUCHOS_ASSERT(dirichletRows.extent(0) == rowMap->getLocalNumElements());
588  TEUCHOS_ASSERT(dirichletCols.extent(0) == colMap->getLocalNumElements());
589  TEUCHOS_ASSERT(dirichletDomain.extent(0) == domMap->getLocalNumElements());
590  RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, /*zeroOut=*/true);
591  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
592  auto myColsToZeroView = myColsToZero->getDeviceLocalView(Xpetra::Access::ReadWrite);
593  auto localMatrix = A.getLocalMatrixDevice();
594  Kokkos::parallel_for("MueLu:Maxwell1::DetectDirichletCols", range_type(0,rowMap->getLocalNumElements()),
595  KOKKOS_LAMBDA(const LocalOrdinal row) {
596  if (dirichletRows(row)) {
597  auto rowView = localMatrix.row(row);
598  auto length = rowView.length;
599 
600  for (decltype(length) colID = 0; colID < length; colID++)
601  myColsToZeroView(rowView.colidx(colID),0) = impl_ATS::one();
602  }
603  });
604 
605  RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
606  RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
607  if (!importer.is_null()) {
608  globalColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, /*zeroOut=*/true);
609  // export to domain map
610  globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
611  // import to column map
612  myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
613  }
614  else
615  globalColsToZero = myColsToZero;
616  FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(globalColsToZero->getDeviceLocalView(Xpetra::Access::ReadOnly),dirichletDomain);
617  FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(myColsToZero->getDeviceLocalView(Xpetra::Access::ReadOnly),dirichletCols);
618  }
619 
620 
621  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
622  void
623  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
624  DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
625  const Kokkos::View<bool*, typename Node::device_type> & dirichletRows,
626  Kokkos::View<bool*, typename Node::device_type> dirichletCols,
627  Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
628  MueLu::DetectDirichletColsAndDomains<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A,dirichletRows,dirichletCols,dirichletDomain);
629  }
630 
631  template <class Node>
632  void
633  Utilities_kokkos<double,int,int,Node>::
634  DetectDirichletColsAndDomains(const Xpetra::Matrix<double,int,int,Node>& A,
635  const Kokkos::View<bool*, typename Node::device_type> & dirichletRows,
636  Kokkos::View<bool*, typename Node::device_type> dirichletCols,
637  Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
638  MueLu::DetectDirichletColsAndDomains<double,int,int,Node>(A,dirichletRows,dirichletCols,dirichletDomain);
639  }
640 
641 
642  // Zeros out rows
643  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
644  void
645  ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& A,
646  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
647  Scalar replaceWith) {
648  using ATS = Kokkos::ArithTraits<Scalar>;
649  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
650 
651  typename ATS::val_type impl_replaceWith = replaceWith;
652 
653  auto localMatrix = A->getLocalMatrixDevice();
654  LocalOrdinal numRows = A->getLocalNumRows();
655 
656  Kokkos::parallel_for("MueLu:Utils::ZeroDirichletRows", range_type(0,numRows),
657  KOKKOS_LAMBDA(const LocalOrdinal row) {
658  if (dirichletRows(row)) {
659  auto rowView = localMatrix.row(row);
660  auto length = rowView.length;
661  for (decltype(length) colID = 0; colID < length; colID++)
662  rowView.value(colID) = impl_replaceWith;
663  }
664  });
665  }
666 
667  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
668  void
669  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
670  ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
671  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows,
672  Scalar replaceWith) {
673  MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows, replaceWith);
674  }
675 
676  template <class Node>
677  void
678  Utilities_kokkos<double,int,int,Node>::
679  ZeroDirichletRows(RCP<Xpetra::Matrix<double, int, int, Node> >& A,
680  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
681  double replaceWith) {
682  return MueLu::ZeroDirichletRows<double,int,int,Node>(A, dirichletRows, replaceWith);
683  }
684 
685 
686  // Zeros out rows
687  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
688  void
689  ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
690  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
691  Scalar replaceWith) {
692  using ATS = Kokkos::ArithTraits<Scalar>;
693  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
694 
695  typename ATS::val_type impl_replaceWith = replaceWith;
696 
697  auto myCols = X->getDeviceLocalView(Xpetra::Access::ReadWrite);
698  size_t numVecs = X->getNumVectors();
699  Kokkos::parallel_for("MueLu:Utils::ZeroDirichletRows_MV", range_type(0,dirichletRows.size()),
700  KOKKOS_LAMBDA(const size_t i) {
701  if (dirichletRows(i)) {
702  for(size_t j=0; j<numVecs; j++)
703  myCols(i,j) = impl_replaceWith;
704  }
705  });
706  }
707 
708  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
709  void
710  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
711  ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >& X,
712  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows,
713  Scalar replaceWith) {
714  MueLu::ZeroDirichletRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X, dirichletRows, replaceWith);
715  }
716 
717  template <class Node>
718  void
719  Utilities_kokkos<double,int,int,Node>::
720  ZeroDirichletRows(RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, int, int, Node> >& X,
721  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
722  double replaceWith) {
723  return MueLu::ZeroDirichletRows<double,int,int,Node>(X, dirichletRows, replaceWith);
724  }
725 
726 
727  // Zeros out columns
728  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
729  void
730  ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
731  const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols,
732  Scalar replaceWith) {
733  using ATS = Kokkos::ArithTraits<Scalar>;
734  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
735 
736  typename ATS::val_type impl_replaceWith = replaceWith;
737 
738  auto localMatrix = A->getLocalMatrixDevice();
739  LocalOrdinal numRows = A->getLocalNumRows();
740 
741  Kokkos::parallel_for("MueLu:Utils::ZeroDirichletCols", range_type(0,numRows),
742  KOKKOS_LAMBDA(const LocalOrdinal row) {
743  auto rowView = localMatrix.row(row);
744  auto length = rowView.length;
745  for (decltype(length) colID = 0; colID < length; colID++)
746  if (dirichletCols(rowView.colidx(colID))) {
747  rowView.value(colID) = impl_replaceWith;
748  }
749  });
750  }
751 
752  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
753  void
754  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
755  ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
756  const Kokkos::View<const bool*, typename NO::device_type>& dirichletCols,
757  Scalar replaceWith) {
758  MueLu::ZeroDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletCols, replaceWith);
759  }
760 
761  template <class Node>
762  void
763  Utilities_kokkos<double,int,int,Node>::
764  ZeroDirichletCols(RCP<Xpetra::Matrix<double,int,int,Node> >& A,
765  const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols,
766  double replaceWith) {
767  return MueLu::ZeroDirichletCols<double,int,int,Node>(A, dirichletCols, replaceWith);
768  }
769 
770  // Applies rowsum criterion
771  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
772  void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
773  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
774  Kokkos::View<bool*, typename Node::device_type> & dirichletRows)
775  {
776  typedef Teuchos::ScalarTraits<Scalar> STS;
777  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> rowmap = A.getRowMap();
778 
779  auto dirichletRowsHost = Kokkos::create_mirror_view(dirichletRows);
780  Kokkos::deep_copy(dirichletRowsHost, dirichletRows);
781 
782  for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
783  size_t nnz = A.getNumEntriesInLocalRow(row);
784  ArrayView<const LocalOrdinal> indices;
785  ArrayView<const Scalar> vals;
786  A.getLocalRowView(row, indices, vals);
787 
788  Scalar rowsum = STS::zero();
789  Scalar diagval = STS::zero();
790  for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
791  LocalOrdinal col = indices[colID];
792  if (row == col)
793  diagval = vals[colID];
794  rowsum += vals[colID];
795  }
796  if (STS::real(rowsum) > STS::magnitude(diagval) * rowSumTol)
797  dirichletRowsHost(row) = true;
798  }
799 
800  Kokkos::deep_copy(dirichletRows, dirichletRowsHost);
801  }
802 
803  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
804  void
805  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
806  ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
807  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
808  Kokkos::View<bool*, typename Node::device_type> & dirichletRows)
809  {
810  MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A,rowSumTol,dirichletRows);
811  }
812 
813 
814  template <class Node>
815  void
816  Utilities_kokkos<double,int,int,Node>::
817  ApplyRowSumCriterion(const Xpetra::Matrix<double,int,int,Node>& A,
818  const typename Teuchos::ScalarTraits<double>::magnitudeType rowSumTol,
819  Kokkos::View<bool*, typename Node::device_type> & dirichletRows)
820  {
821  MueLu::ApplyRowSumCriterion<double, int, int, Node>(A,rowSumTol,dirichletRows);
822  }
823 
824 
825  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
826  RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
827  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
828  RealValuedToScalarMultiVector(RCP<RealValuedMultiVector > X) {
829  RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xscalar;
830 #if defined(HAVE_XPETRA_TPETRA) && (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))
831  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
832 
833  // Need to cast the real-valued multivector to Scalar=complex
834  if ((typeid(Scalar).name() == typeid(std::complex<double>).name()) ||
835  (typeid(Scalar).name() == typeid(std::complex<float>).name())) {
836  size_t numVecs = X->getNumVectors();
837  Xscalar = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X->getMap(),numVecs);
838  auto XVec = X->getDeviceLocalView(Xpetra::Access::ReadOnly);
839  auto XVecScalar = Xscalar->getDeviceLocalView(Xpetra::Access::ReadWrite);
840 
841  Kokkos::parallel_for("MueLu:Utils::RealValuedToScalarMultiVector", range_type(0,X->getLocalLength()),
842  KOKKOS_LAMBDA(const size_t i) {
843  for (size_t j=0; j<numVecs; j++)
844  XVecScalar(i,j) = XVec(i,j);
845  });
846  } else
847 #endif
848  Xscalar = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(X);
849  return Xscalar;
850  }
851 
852  template <class Node>
853  RCP<Xpetra::MultiVector<double,int,int,Node> >
854  Utilities_kokkos<double,int,int,Node>::
855  RealValuedToScalarMultiVector(RCP<Xpetra::MultiVector<Magnitude,int,int,Node> > X) {
856  return X;
857  }
858 
859 
860  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
861  Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > ReverseCuthillMcKee(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Op) {
862  using local_matrix_type = typename Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::local_matrix_type;
863  using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
864  using lno_nnz_view_t = typename local_graph_type::entries_type::non_const_type;
865  using device = typename local_graph_type::device_type;
866  using execution_space = typename local_matrix_type::execution_space;
867  using ordinal_type = typename local_matrix_type::ordinal_type;
868 
869  local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
870 
871  lno_nnz_view_t rcmOrder = KokkosGraph::Experimental::graph_rcm
872  <device, typename local_graph_type::row_map_type, typename local_graph_type::entries_type, lno_nnz_view_t>
873  (localGraph.row_map, localGraph.entries);
874 
875  RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > retval =
876  Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRowMap());
877 
878  // Copy out and reorder data
879  auto view1D = Kokkos::subview(retval->getDeviceLocalView(Xpetra::Access::ReadWrite),Kokkos::ALL (), 0);
880  Kokkos::parallel_for("Utilities_kokkos::ReverseCuthillMcKee",
881  Kokkos::RangePolicy<ordinal_type, execution_space>(0, localGraph.numRows()),
882  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
883  view1D(rcmOrder(rowIdx)) = rowIdx;
884  });
885  return retval;
886  }
887 
888  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
889  Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > CuthillMcKee(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Op) {
890  using local_matrix_type = typename Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::local_matrix_type;
891  using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
892  using lno_nnz_view_t = typename local_graph_type::entries_type::non_const_type;
893  using device = typename local_graph_type::device_type;
894  using execution_space = typename local_matrix_type::execution_space;
895  using ordinal_type = typename local_matrix_type::ordinal_type;
896 
897  local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
898  LocalOrdinal numRows = localGraph.numRows();
899 
900  lno_nnz_view_t rcmOrder = KokkosGraph::Experimental::graph_rcm
901  <device, typename local_graph_type::row_map_type, typename local_graph_type::entries_type, lno_nnz_view_t>
902  (localGraph.row_map, localGraph.entries);
903 
904  RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > retval =
905  Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRowMap());
906 
907  // Copy out data
908  auto view1D = Kokkos::subview(retval->getDeviceLocalView(Xpetra::Access::ReadWrite),Kokkos::ALL (), 0);
909  // Since KokkosKernels produced RCM, also reverse the order of the view to get CM
910  Kokkos::parallel_for("Utilities_kokkos::ReverseCuthillMcKee",
911  Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
912  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
913  view1D(rcmOrder(numRows - 1 - rowIdx)) = rowIdx;
914  });
915  return retval;
916  }
917 
918  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
919  Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> >
920  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReverseCuthillMcKee(const Matrix &Op) {
921  return MueLu::ReverseCuthillMcKee<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Op);
922  }
923 
924  template <class Node>
925  Teuchos::RCP<Xpetra::Vector<int,int,int,Node> >
926  Utilities_kokkos<double,int,int,Node>::ReverseCuthillMcKee(const Matrix &Op) {
927  return MueLu::ReverseCuthillMcKee<double,int,int,Node>(Op);
928  }
929 
930  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
931  Teuchos::RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> >
932  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CuthillMcKee(const Matrix &Op) {
933  return MueLu::CuthillMcKee<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Op);
934  }
935 
936  template <class Node>
937  Teuchos::RCP<Xpetra::Vector<int,int,int,Node> >
938  Utilities_kokkos<double,int,int,Node>::CuthillMcKee(const Matrix &Op) {
939  return MueLu::CuthillMcKee<double,int,int,Node>(Op);
940  }
941 
942  // Applies Ones-and-Zeros to matrix rows
943  // Takes a Boolean array.
944  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
945  void
946  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
947  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
948  TEUCHOS_ASSERT(A->isFillComplete());
949  using ATS = Kokkos::ArithTraits<Scalar>;
950  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
951  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
952 
953  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A->getDomainMap();
954  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > ranMap = A->getRangeMap();
955  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Rmap = A->getRowMap();
956  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Cmap = A->getColMap();
957 
958  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getLocalNumElements());
959 
960  auto localMatrix = A->getLocalMatrixDevice();
961  auto localRmap = Rmap->getLocalMap();
962  auto localCmap = Cmap->getLocalMap();
963 
964  Kokkos::parallel_for("MueLu::Utils::ApplyOAZ",range_type(0,dirichletRows.extent(0)),
965  KOKKOS_LAMBDA(const LocalOrdinal row) {
966  if (dirichletRows(row)){
967  auto rowView = localMatrix.row(row);
968  auto length = rowView.length;
969  auto row_gid = localRmap.getGlobalElement(row);
970  auto row_lid = localCmap.getLocalElement(row_gid);
971 
972  for (decltype(length) colID = 0; colID < length; colID++)
973  if (rowView.colidx(colID) == row_lid)
974  rowView.value(colID) = impl_ATS::one();
975  else
976  rowView.value(colID) = impl_ATS::zero();
977  }
978  });
979  }
980 
981  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
982  void
983  Utilities_kokkos<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
984  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
985  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
986  MueLu::ApplyOAZToMatrixRows<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, dirichletRows);
987  }
988 
989  template <class Node>
990  void
991  Utilities_kokkos<double,int,int,Node>::
992  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<double,int,int,Node> >& A,
993  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
994  MueLu::ApplyOAZToMatrixRows<double,int,int,Node>(A, dirichletRows);
995  }
996 
997 } //namespace MueLu
998 
999 #define MUELU_UTILITIES_KOKKOS_SHORT
1000 #endif // MUELU_UTILITIES_KOKKOS_DEF_HPP
Xpetra::ExportFactory::Build
static RCP< Export< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
Xpetra_Vector.hpp
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:50
MueLu::Utilities_kokkos::DetectDirichletCols
static Kokkos::View< bool *, typename NO::device_type > DetectDirichletCols(const Matrix &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
Definition: MueLu_Utilities_kokkos_def.hpp:529
MueLu::Utilities_kokkos::MyOldScaleMatrix
static void MyOldScaleMatrix(Matrix &Op, const Teuchos::ArrayRCP< const SC > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
Definition: MueLu_Utilities_kokkos_def.hpp:250
EpetraExt_CrsMatrixIn.h
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
MueLu::ReverseCuthillMcKee
Teuchos::RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Definition: MueLu_Utilities_kokkos_def.hpp:861
Tpetra::CrsMatrix::getRangeMap
Teuchos::RCP< const map_type > getRangeMap() const override
GO
GlobalOrdinal GO
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Tpetra::CrsMatrix::resumeFill
void resumeFill(const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
MueLu::ZeroDirichletCols
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
Definition: MueLu_Utilities_kokkos_def.hpp:730
Xpetra_EpetraUtils.hpp
MueLu::ZeroDirichletRows
void ZeroDirichletRows(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, Scalar replaceWith)
Definition: MueLu_Utilities_kokkos_def.hpp:645
LO
LocalOrdinal LO
MueLu::CuthillMcKee
Teuchos::RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Definition: MueLu_Utilities_kokkos_def.hpp:889
Tpetra::BlockCrsMatrix::getLocalNumRows
size_t getLocalNumRows() const override
MueLu::Utilities_kokkos::RealValuedToScalarMultiVector
static RCP< MultiVector > RealValuedToScalarMultiVector(RCP< RealValuedMultiVector > X)
Definition: MueLu_Utilities_kokkos_def.hpp:828
MueLu::ApplyOAZToMatrixRows
void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows)
Definition: MueLu_Utilities_kokkos_def.hpp:946
EpetraExt_MatrixMatrix.h
Xpetra_Matrix.hpp
Teuchos::ArrayRCP::size
size_type size() const
Xpetra_MatrixMatrix.hpp
MueLu_ConfigDefs.hpp
MueLu_Utilities_kokkos_decl.hpp
MueLu::Utilities_kokkos::GetMatrixDiagonal
static RCP< Vector > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
Definition: MueLu_Utilities_kokkos_def.hpp:117
Teuchos::TimeMonitor::getNewTimer
static RCP< Time > getNewTimer(const std::string &name)
MueLu::Utilities_kokkos::GetMatrixDiagonalInverse
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=TST::eps()*100, const bool doLumped=false)
Extract Matrix Diagonal.
Definition: MueLu_Utilities_kokkos_def.hpp:204
MueLu::DetectDirichletColsAndDomains
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)
Definition: MueLu_Utilities_kokkos_def.hpp:577
EpetraExt_RowMatrixOut.h
Xpetra::Matrix::fillComplete
virtual void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > &params=null)=0
Xpetra::Matrix::getLocalDiagCopy
virtual void getLocalDiagCopy(Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &diag) const =0
Tpetra::CrsMatrix::getLocalMaxNumRowEntries
size_t getLocalMaxNumRowEntries() const override
Xpetra::CrsMatrixWrap::getLocalDiagOffsets
void getLocalDiagOffsets(Teuchos::ArrayRCP< size_t > &offsets) const
Xpetra::Matrix::getLocalRowView
virtual void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const =0
Xpetra_MultiVector.hpp
MueLu::DetectDirichletRows
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)
Definition: MueLu_Utilities_kokkos_def.hpp:365
Tpetra::CrsMatrix::isFillComplete
bool isFillComplete() const override
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos_ParameterList.hpp
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:49
EpetraExt_BlockMapOut.h
Xpetra::CrsMatrixWrap::getLocalDiagCopy
void getLocalDiagCopy(Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &diag) const
EpetraExt_MultiVectorOut.h
MueLu::Utilities_kokkos::GetMatrixOverlappedDiagonal
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
Definition: MueLu_Utilities_kokkos_def.hpp:220
EpetraExt_MultiVectorIn.h
Tpetra::CrsMatrix::replaceGlobalValues
local_ordinal_type replaceGlobalValues(const global_ordinal_type globalRow, const Kokkos::View< const global_ordinal_type *, Kokkos::AnonymousSpace > &inputInds, const Kokkos::View< const impl_scalar_type *, Kokkos::AnonymousSpace > &inputVals)
MueLu::Utilities_kokkos::ApplyRowSumCriterion
static void ApplyRowSumCriterion(const Matrix &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, typename NO::device_type > &dirichletRows)
Definition: MueLu_Utilities_kokkos_def.hpp:806
Tpetra::CrsMatrix::getNumEntriesInLocalRow
size_t getNumEntriesInLocalRow(local_ordinal_type localRow) const override
Xpetra::Matrix::getNumEntriesInLocalRow
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
EpetraExt_BlockMapIn.h
MueLu::Utilities_kokkos< double, int, int, Node >::Magnitude
Teuchos::ScalarTraits< Scalar >::magnitudeType Magnitude
Definition: MueLu_Utilities_kokkos_decl.hpp:356
MueLu::GetMatrixDiagonalInverse
Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixDiagonalInverse(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, typename Teuchos::ScalarTraits< Scalar >::magnitudeType tol, const bool doLumped)
Definition: MueLu_Utilities_kokkos_def.hpp:128
MueLu::Utilities_kokkos::DetectDirichletColsAndDomains
static 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)
Detects Dirichlet columns &amp; domains from a list of Dirichlet rows.
Definition: MueLu_Utilities_kokkos_def.hpp:624
Xpetra::VectorFactory::Build
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
MueLu::Utilities_kokkos::FindNonZeros
static 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)
Find non-zero values in an ArrayRCP Compares the value to 2 * machine epsilon.
Definition: MueLu_Utilities_kokkos_def.hpp:562
MueLu::Utilities_kokkos::Magnitude
typename TST::magnitudeType Magnitude
Definition: MueLu_Utilities_kokkos_decl.hpp:114
Xpetra_EpetraMap.hpp
Xpetra_BlockedCrsMatrix.hpp
Xpetra_MatrixFactory.hpp
MueLu::Utilities_kokkos::DetectDirichletRows
static Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
Definition: MueLu_Utilities_kokkos_def.hpp:463
Xpetra_DefaultPlatform.hpp
Xpetra::Matrix::getCrsGraph
virtual RCP< const CrsGraph > getCrsGraph() const =0
Tpetra::CrsMatrix::getDomainMap
Teuchos::RCP< const map_type > getDomainMap() const override
Xpetra::Matrix::isLocallyIndexed
virtual bool isLocallyIndexed() const =0
Xpetra_EpetraMultiVector.hpp
MueLu::FindNonZeros
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)
Definition: MueLu_Utilities_kokkos_def.hpp:545
SC
Scalar SC
Tpetra::BlockCrsMatrix::getBlockSize
LO getBlockSize() const
Tpetra::CrsMatrix::getGlobalRowView
virtual void getGlobalRowView(GlobalOrdinal GlobalRow, global_inds_host_view_type &indices, values_host_view_type &values) const =0
MueLu::ApplyRowSumCriterion
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)
Definition: MueLu_Utilities_kokkos_def.hpp:772
Tpetra::BlockCrsMatrix::getValuesDevice
impl_scalar_type_dualview::t_dev::const_type getValuesDevice(const LO &lclRow) const
Xpetra::Operator::getRangeMap
virtual Teuchos::RCP< const Map > getRangeMap() const =0
MueLu::Utilities_kokkos::ReverseCuthillMcKee
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Matrix &Op)
Definition: MueLu_Utilities_kokkos_def.hpp:920
Xpetra_VectorFactory.hpp
Tpetra::CrsMatrix::getNumEntriesInGlobalRow
size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const override
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MatrixMarket_Tpetra.hpp
Tpetra::CrsMatrix::getLocalRowView
virtual void getLocalRowView(LocalOrdinal LocalRow, local_inds_host_view_type &indices, values_host_view_type &values) const =0
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Tpetra::CrsMatrix::replaceLocalValues
local_ordinal_type replaceLocalValues(const local_ordinal_type localRow, const Kokkos::View< const local_ordinal_type *, Kokkos::AnonymousSpace > &inputInds, const Kokkos::View< const impl_scalar_type *, Kokkos::AnonymousSpace > &inputVals)
MueLu::Utilities_kokkos::MyOldScaleMatrix_Epetra
static void MyOldScaleMatrix_Epetra(Matrix &Op, const Teuchos::ArrayRCP< SC > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
Definition: MueLu_Utilities_kokkos_def.hpp:285
MueLu::Utilities_kokkos::ApplyOAZToMatrixRows
static void ApplyOAZToMatrixRows(RCP< Matrix > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows)
Definition: MueLu_Utilities_kokkos_def.hpp:984
MueLu::DetectDirichletCols
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)
Definition: MueLu_Utilities_kokkos_def.hpp:477
Xpetra_Import.hpp
MueLu::Utilities_kokkos::MyOldScaleMatrix_Tpetra
static void MyOldScaleMatrix_Tpetra(Matrix &Op, const Teuchos::ArrayRCP< SC > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
Definition: MueLu_Utilities_kokkos_def.hpp:290
Xpetra_Operator.hpp
MueLu::Utilities_kokkos::ZeroDirichletCols
static void ZeroDirichletCols(RCP< Matrix > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletCols, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
Definition: MueLu_Utilities_kokkos_def.hpp:755
Xpetra::Operator::getDomainMap
virtual Teuchos::RCP< const Map > getDomainMap() const =0
MueLu::Utilities_kokkos::ZeroDirichletRows
static void ZeroDirichletRows(RCP< Matrix > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
Xpetra_ImportFactory.hpp
MueLu::Utilities_kokkos::CuthillMcKee
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Matrix &Op)
Definition: MueLu_Utilities_kokkos_def.hpp:932
Teuchos::RCP::is_null
bool is_null() const
Tpetra::CrsMatrix::getRowMap
Teuchos::RCP< const map_type > getRowMap() const override