MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_UtilitiesBase_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_UTILITIESBASE_DEF_HPP
11 #define MUELU_UTILITIESBASE_DEF_HPP
12 
13 #include "MueLu_ConfigDefs.hpp"
14 
15 #include "MueLu_UtilitiesBase_decl.hpp"
16 
17 #include <Kokkos_Core.hpp>
18 #include <KokkosSparse_CrsMatrix.hpp>
19 #include <KokkosSparse_getDiagCopy.hpp>
20 
21 #include <Xpetra_BlockedVector.hpp>
22 #include <Xpetra_BlockedMap.hpp>
23 #include <Xpetra_BlockedMultiVector.hpp>
24 #include <Xpetra_ExportFactory.hpp>
25 
26 #include <Xpetra_Import.hpp>
27 #include <Xpetra_ImportFactory.hpp>
28 #include <Xpetra_MatrixMatrix.hpp>
29 #include <Xpetra_CrsGraph.hpp>
30 #include <Xpetra_CrsGraphFactory.hpp>
31 #include <Xpetra_CrsMatrixWrap.hpp>
32 #include <Xpetra_StridedMap.hpp>
33 
34 #include "MueLu_Exceptions.hpp"
35 #include "Xpetra_CrsMatrixFactory.hpp"
36 
37 #include <KokkosKernels_Handle.hpp>
38 #include <KokkosGraph_RCM.hpp>
39 
40 namespace MueLu {
41 
42 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
43 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
44 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
45  Crs2Op(RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> Op) {
46  if (Op.is_null())
47  return Teuchos::null;
48  return rcp(new CrsMatrixWrap(Op));
49 }
50 
51 template <class Scalar,
52  class LocalOrdinal,
53  class GlobalOrdinal,
54  class Node>
55 Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
56 removeSmallEntries(Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
57  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType threshold,
58  const bool keepDiagonal) {
59  using crs_matrix = Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
60  using row_ptr_type = typename crs_matrix::local_graph_type::row_map_type::non_const_type;
61  using col_idx_type = typename crs_matrix::local_graph_type::entries_type::non_const_type;
62  using vals_type = typename crs_matrix::local_matrix_type::values_type;
63  using execution_space = typename crs_matrix::local_matrix_type::execution_space;
64 
65 #if KOKKOS_VERSION >= 40799
66  using ATS = KokkosKernels::ArithTraits<Scalar>;
67 #else
68  using ATS = Kokkos::ArithTraits<Scalar>;
69 #endif
70  using impl_SC = typename ATS::val_type;
71 #if KOKKOS_VERSION >= 40799
72  using impl_ATS = KokkosKernels::ArithTraits<impl_SC>;
73 #else
74  using impl_ATS = Kokkos::ArithTraits<impl_SC>;
75 #endif
76 
77  auto lclA = A->getLocalMatrixDevice();
78 
79  auto rowptr = row_ptr_type("rowptr", lclA.numRows() + 1);
80  col_idx_type idx;
81  vals_type vals;
82  LocalOrdinal nnz;
83 
84  if (keepDiagonal) {
85  auto lclRowMap = A->getRowMap()->getLocalMap();
86  auto lclColMap = A->getColMap()->getLocalMap();
87  Kokkos::parallel_scan(
88  "removeSmallEntries::rowptr",
89  Kokkos::RangePolicy<LocalOrdinal, execution_space>(0, lclA.numRows()),
90  KOKKOS_LAMBDA(const LocalOrdinal rlid, LocalOrdinal& partial_nnz, bool is_final) {
91  auto row = lclA.row(rlid);
92  auto rowInCol = lclColMap.getLocalElement(lclRowMap.getGlobalElement(rlid));
93  for (LocalOrdinal k = 0; k < row.length; ++k) {
94  if ((impl_ATS::magnitude(row.value(k)) > threshold) || (row.colidx(k) == rowInCol)) {
95  partial_nnz += 1;
96  }
97  }
98  if (is_final)
99  rowptr(rlid + 1) = partial_nnz;
100  },
101  nnz);
102 
103  idx = col_idx_type("idx", nnz);
104  vals = vals_type("vals", nnz);
105 
106  Kokkos::parallel_for(
107  "removeSmallEntries::indicesValues",
108  Kokkos::RangePolicy<LocalOrdinal, execution_space>(0, lclA.numRows()),
109  KOKKOS_LAMBDA(const LocalOrdinal rlid) {
110  auto row = lclA.row(rlid);
111  auto rowInCol = lclColMap.getLocalElement(lclRowMap.getGlobalElement(rlid));
112  auto I = rowptr(rlid);
113  for (LocalOrdinal k = 0; k < row.length; ++k) {
114  if ((impl_ATS::magnitude(row.value(k)) > threshold) || (row.colidx(k) == rowInCol)) {
115  idx(I) = row.colidx(k);
116  vals(I) = row.value(k);
117  I += 1;
118  }
119  }
120  });
121 
122  Kokkos::fence();
123  } else {
124  Kokkos::parallel_scan(
125  "removeSmallEntries::rowptr",
126  Kokkos::RangePolicy<LocalOrdinal, execution_space>(0, lclA.numRows()),
127  KOKKOS_LAMBDA(const LocalOrdinal rlid, LocalOrdinal& partial_nnz, bool is_final) {
128  auto row = lclA.row(rlid);
129  for (LocalOrdinal k = 0; k < row.length; ++k) {
130  if (impl_ATS::magnitude(row.value(k)) > threshold) {
131  partial_nnz += 1;
132  }
133  }
134  if (is_final)
135  rowptr(rlid + 1) = partial_nnz;
136  },
137  nnz);
138 
139  idx = col_idx_type("idx", nnz);
140  vals = vals_type("vals", nnz);
141 
142  Kokkos::parallel_for(
143  "removeSmallEntries::indicesValues",
144  Kokkos::RangePolicy<LocalOrdinal, execution_space>(0, lclA.numRows()),
145  KOKKOS_LAMBDA(const LocalOrdinal rlid) {
146  auto row = lclA.row(rlid);
147  auto I = rowptr(rlid);
148  for (LocalOrdinal k = 0; k < row.length; ++k) {
149  if (impl_ATS::magnitude(row.value(k)) > threshold) {
150  idx(I) = row.colidx(k);
151  vals(I) = row.value(k);
152  I += 1;
153  }
154  }
155  });
156 
157  Kokkos::fence();
158  }
159 
160  auto lclNewA = typename crs_matrix::local_matrix_type("thresholdedMatrix", lclA.numRows(), lclA.numCols(), nnz, vals, rowptr, idx);
161  auto newA = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(lclNewA, A->getRowMap(), A->getColMap(), A->getDomainMap(), A->getRangeMap());
162 
163  return newA;
164 }
165 
166 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
167 RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
168 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
169  GetThresholdedMatrix(const RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& Ain, const typename Teuchos::ScalarTraits<Scalar>::magnitudeType threshold, const bool keepDiagonal, const GlobalOrdinal expectedNNZperRow) {
170  auto crsWrap = rcp_dynamic_cast<CrsMatrixWrap>(Ain);
171  if (!crsWrap.is_null()) {
172  auto crsMat = crsWrap->getCrsMatrix();
173  auto filteredMat = removeSmallEntries(crsMat, threshold, keepDiagonal);
174  return rcp_static_cast<CrsMatrixWrap>(filteredMat);
175  }
176 
177  RCP<const Map> rowmap = Ain->getRowMap();
178  RCP<const Map> colmap = Ain->getColMap();
179  RCP<CrsMatrixWrap> Aout = rcp(new CrsMatrixWrap(rowmap, expectedNNZperRow <= 0 ? Ain->getGlobalMaxNumRowEntries() : expectedNNZperRow));
180  // loop over local rows
181  for (size_t row = 0; row < Ain->getLocalNumRows(); row++) {
182  size_t nnz = Ain->getNumEntriesInLocalRow(row);
183 
184  Teuchos::ArrayView<const LocalOrdinal> indices;
185  Teuchos::ArrayView<const Scalar> vals;
186  Ain->getLocalRowView(row, indices, vals);
187 
188  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz, Exceptions::RuntimeError, "MueLu::ThresholdAFilterFactory::Build: number of nonzeros not equal to number of indices? Error.");
189 
190  Teuchos::ArrayRCP<GlobalOrdinal> indout(indices.size(), Teuchos::ScalarTraits<GlobalOrdinal>::zero());
191  Teuchos::ArrayRCP<Scalar> valout(indices.size(), Teuchos::ScalarTraits<Scalar>::zero());
192  size_t nNonzeros = 0;
193  if (keepDiagonal) {
194  GlobalOrdinal glbRow = rowmap->getGlobalElement(row);
195  LocalOrdinal lclColIdx = colmap->getLocalElement(glbRow);
196  for (size_t i = 0; i < (size_t)indices.size(); i++) {
197  if (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(threshold) || indices[i] == lclColIdx) {
198  indout[nNonzeros] = colmap->getGlobalElement(indices[i]); // LID -> GID (column)
199  valout[nNonzeros] = vals[i];
200  nNonzeros++;
201  }
202  }
203  } else
204  for (size_t i = 0; i < (size_t)indices.size(); i++) {
205  if (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(threshold)) {
206  indout[nNonzeros] = colmap->getGlobalElement(indices[i]); // LID -> GID (column)
207  valout[nNonzeros] = vals[i];
208  nNonzeros++;
209  }
210  }
211 
212  indout.resize(nNonzeros);
213  valout.resize(nNonzeros);
214 
215  Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0, indout.size()), valout.view(0, valout.size()));
216  }
217  Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());
218 
219  return Aout;
220 }
221 
222 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
223 RCP<Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node>>
224 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
225  GetThresholdedGraph(const RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow) {
226  using STS = Teuchos::ScalarTraits<Scalar>;
227  RCP<CrsGraph> sparsityPattern = CrsGraphFactory::Build(A->getRowMap(), expectedNNZperRow <= 0 ? A->getGlobalMaxNumRowEntries() : expectedNNZperRow);
228 
229  RCP<Vector> diag = GetMatrixOverlappedDiagonal(*A);
230  ArrayRCP<const Scalar> D = diag->getData(0);
231 
232  for (size_t row = 0; row < A->getLocalNumRows(); row++) {
233  ArrayView<const LocalOrdinal> indices;
234  ArrayView<const Scalar> vals;
235  A->getLocalRowView(row, indices, vals);
236 
237  GlobalOrdinal globalRow = A->getRowMap()->getGlobalElement(row);
238  LocalOrdinal col = A->getColMap()->getLocalElement(globalRow);
239 
240  const Scalar Dk = STS::magnitude(D[col]) > 0.0 ? STS::magnitude(D[col]) : 1.0;
241  Array<GlobalOrdinal> indicesNew;
242 
243  for (size_t i = 0; i < size_t(indices.size()); i++)
244  // keep diagonal per default
245  if (col == indices[i] || STS::magnitude(STS::squareroot(Dk) * vals[i] * STS::squareroot(Dk)) > STS::magnitude(threshold))
246  indicesNew.append(A->getColMap()->getGlobalElement(indices[i]));
247 
248  sparsityPattern->insertGlobalIndices(globalRow, ArrayView<const GlobalOrdinal>(indicesNew.data(), indicesNew.length()));
249  }
250  sparsityPattern->fillComplete();
251 
252  return sparsityPattern;
253 }
254 
255 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
256 Teuchos::ArrayRCP<Scalar>
257 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
258  GetMatrixDiagonal_arcp(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A) {
259  size_t numRows = A.getRowMap()->getLocalNumElements();
260  Teuchos::ArrayRCP<Scalar> diag(numRows);
261  Teuchos::ArrayView<const LocalOrdinal> cols;
262  Teuchos::ArrayView<const Scalar> vals;
263  for (size_t i = 0; i < numRows; ++i) {
264  A.getLocalRowView(i, cols, vals);
265  LocalOrdinal j = 0;
266  for (; j < cols.size(); ++j) {
267  if (Teuchos::as<size_t>(cols[j]) == i) {
268  diag[i] = vals[j];
269  break;
270  }
271  }
272  if (j == cols.size()) {
273  // Diagonal entry is absent
274  diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
275  }
276  }
277  return diag;
278 }
279 
280 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
281 Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
282 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
283  GetMatrixDiagonal(const Matrix& A) {
284  const auto rowMap = A.getRowMap();
285  auto diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, true);
286 
287  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
288  if ((crsOp != NULL) && (rowMap->lib() == Xpetra::UseTpetra)) {
289  using device_type = typename CrsGraph::device_type;
290  Kokkos::View<size_t*, device_type> offsets("offsets", rowMap->getLocalNumElements());
291  crsOp->getCrsGraph()->getLocalDiagOffsets(offsets);
292  crsOp->getCrsMatrix()->getLocalDiagCopy(*diag, offsets);
293  } else {
294  A.getLocalDiagCopy(*diag);
295  }
296 
297  return diag;
298 }
299 
300 // template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
301 // RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
302 // UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
303 // GetMatrixDiagonalInverse(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> & A, Magnitude tol, Scalar valReplacement) {
304 // Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer("UtilitiesBase::GetMatrixDiagonalInverse");
305 
306 // RCP<const Map> rowMap = A.getRowMap();
307 // RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
308 
309 // A.getLocalDiagCopy(*diag);
310 
311 // RCP<Vector> inv = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(diag, tol, valReplacement);
312 
313 // return inv;
314 // }
315 
316 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
317 Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
318 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
319  GetMatrixDiagonalInverse(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
320  typename Teuchos::ScalarTraits<Scalar>::magnitudeType tol,
321  Scalar valReplacement,
322  const bool doLumped) {
323  Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer("Utilities::GetMatrixDiagonalInverse");
324 
325  RCP<const BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(rcpFromRef(A));
326  if (!bA.is_null()) {
327  RCP<const Map> rowMap = A.getRowMap();
328  RCP<Vector> diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, true);
329  A.getLocalDiagCopy(*diag);
330  RCP<Vector> inv = MueLu::UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetInverse(diag, tol, valReplacement);
331  return inv;
332  }
333 
334  // Some useful type definitions
335  using local_matrix_type = typename Matrix::local_matrix_type;
336  // using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
337  using value_type = typename local_matrix_type::value_type;
338  using values_type = typename local_matrix_type::values_type;
339  using scalar_type = typename values_type::non_const_value_type;
340  using ordinal_type = typename local_matrix_type::ordinal_type;
341  using execution_space = typename local_matrix_type::execution_space;
342  // using memory_space = typename local_matrix_type::memory_space;
343  // Be careful with this one, if using KokkosKernels::ArithTraits<Scalar>
344  // you are likely to run into errors when handling std::complex<>
345  // a good way to work around that is to use the following:
346  // using KAT = KokkosKernels::ArithTraits<KokkosKernels::ArithTraits<Scalar>::val_type> >
347  // here we have: value_type = KokkosKernels::ArithTraits<Scalar>::val_type
348 #if KOKKOS_VERSION >= 40799
349  using KAT = KokkosKernels::ArithTraits<value_type>;
350 #else
351  using KAT = Kokkos::ArithTraits<value_type>;
352 #endif
353 
354  // Get/Create distributed objects
355  RCP<const Map> rowMap = A.getRowMap();
356  RCP<Vector> diag = VectorFactory::Build(rowMap, false);
357 
358  // Now generate local objects
359  local_matrix_type localMatrix = A.getLocalMatrixDevice();
360  auto diagVals = diag->getLocalViewDevice(Xpetra::Access::ReadWrite);
361 
362  ordinal_type numRows = localMatrix.graph.numRows();
363 
364  scalar_type valReplacement_dev = valReplacement;
365 
366  // Note: 2019-11-21, LBV
367  // This could be implemented with a TeamPolicy over the rows
368  // and a TeamVectorRange over the entries in a row if performance
369  // becomes more important here.
370  if (!doLumped)
371  Kokkos::parallel_for(
372  "Utilities::GetMatrixDiagonalInverse",
373  Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
374  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
375  bool foundDiagEntry = false;
376  auto myRow = localMatrix.rowConst(rowIdx);
377  for (ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
378  if (myRow.colidx(entryIdx) == rowIdx) {
379  foundDiagEntry = true;
380  if (KAT::magnitude(myRow.value(entryIdx)) > KAT::magnitude(tol)) {
381  diagVals(rowIdx, 0) = KAT::one() / myRow.value(entryIdx);
382  } else {
383  diagVals(rowIdx, 0) = valReplacement_dev;
384  }
385  break;
386  }
387  }
388 
389  if (!foundDiagEntry) {
390  diagVals(rowIdx, 0) = KAT::zero();
391  }
392  });
393  else
394  Kokkos::parallel_for(
395  "Utilities::GetMatrixDiagonalInverse",
396  Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
397  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
398  auto myRow = localMatrix.rowConst(rowIdx);
399  for (ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
400  diagVals(rowIdx, 0) += KAT::magnitude(myRow.value(entryIdx));
401  }
402  if (KAT::magnitude(diagVals(rowIdx, 0)) > KAT::magnitude(tol))
403  diagVals(rowIdx, 0) = KAT::one() / diagVals(rowIdx, 0);
404  else
405  diagVals(rowIdx, 0) = valReplacement_dev;
406  });
407 
408  return diag;
409 }
410 
411 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
412 Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
413 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
414  GetLumpedMatrixDiagonal(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> const& A, const bool doReciprocal,
415  Magnitude tol,
416  Scalar valReplacement,
417  const bool replaceSingleEntryRowWithZero,
418  const bool useAverageAbsDiagVal) {
419  typedef Teuchos::ScalarTraits<Scalar> TST;
420 
421  RCP<Vector> diag = Teuchos::null;
422  const Scalar zero = TST::zero();
423  const Scalar one = TST::one();
424  const Scalar two = one + one;
425 
426  Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
427 
428  RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bA =
429  Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(rcpA);
430  if (bA == Teuchos::null) {
431  RCP<const Map> rowMap = rcpA->getRowMap();
432  diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, true);
433 
434  if (rowMap->lib() == Xpetra::UnderlyingLib::UseTpetra) {
435  Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer("UtilitiesBase::GetLumpedMatrixDiagonal (Kokkos implementation)");
436  // Implement using Kokkos
437  using local_vector_type = typename Vector::dual_view_type::t_dev_um;
438  using local_matrix_type = typename Matrix::local_matrix_type;
439  using execution_space = typename local_vector_type::execution_space;
440  // using rowmap_type = typename local_matrix_type::row_map_type;
441  // using entries_type = typename local_matrix_type::index_type;
442  using values_type = typename local_matrix_type::values_type;
443  using scalar_type = typename values_type::non_const_value_type;
444 #if KOKKOS_VERSION >= 40799
445  using mag_type = typename KokkosKernels::ArithTraits<scalar_type>::mag_type;
446 #else
447  using mag_type = typename Kokkos::ArithTraits<scalar_type>::mag_type;
448 #endif
449 #if KOKKOS_VERSION >= 40799
450  using KAT_S = typename KokkosKernels::ArithTraits<scalar_type>;
451 #else
452  using KAT_S = typename Kokkos::ArithTraits<scalar_type>;
453 #endif
454 #if KOKKOS_VERSION >= 40799
455  using KAT_M = typename KokkosKernels::ArithTraits<mag_type>;
456 #else
457  using KAT_M = typename Kokkos::ArithTraits<mag_type>;
458 #endif
459  using size_type = typename local_matrix_type::non_const_size_type;
460 
461  local_vector_type diag_dev = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
462  local_matrix_type local_mat_dev = rcpA->getLocalMatrixDevice();
463  Kokkos::RangePolicy<execution_space, int> my_policy(0, static_cast<int>(diag_dev.extent(0)));
464  scalar_type valReplacement_dev = valReplacement;
465 
466  if (doReciprocal) {
467  Kokkos::View<int*, execution_space> nnzPerRow("nnz per rows", diag_dev.extent(0));
468  Kokkos::View<scalar_type*, execution_space> regSum("regSum", diag_dev.extent(0));
469  Kokkos::View<mag_type, execution_space> avgAbsDiagVal_dev("avgAbsDiagVal");
470  Kokkos::View<int, execution_space> numDiagsEqualToOne_dev("numDiagsEqualToOne");
471 
472  {
473  Teuchos::TimeMonitor MMM = *Teuchos::TimeMonitor::getNewTimer("GetLumpedMatrixDiagonal: parallel_for (doReciprocal)");
474  Kokkos::parallel_for(
475  "GetLumpedMatrixDiagonal", my_policy,
476  KOKKOS_LAMBDA(const int rowIdx) {
477  diag_dev(rowIdx, 0) = KAT_S::zero();
478  for (size_type entryIdx = local_mat_dev.graph.row_map(rowIdx);
479  entryIdx < local_mat_dev.graph.row_map(rowIdx + 1);
480  ++entryIdx) {
481  regSum(rowIdx) += local_mat_dev.values(entryIdx);
482  if (KAT_M::zero() < KAT_S::abs(local_mat_dev.values(entryIdx))) {
483  ++nnzPerRow(rowIdx);
484  }
485  diag_dev(rowIdx, 0) += KAT_S::abs(local_mat_dev.values(entryIdx));
486  if (rowIdx == local_mat_dev.graph.entries(entryIdx)) {
487  Kokkos::atomic_add(&avgAbsDiagVal_dev(), KAT_S::abs(local_mat_dev.values(entryIdx)));
488  }
489  }
490 
491  if (nnzPerRow(rowIdx) == 1 && KAT_S::magnitude(diag_dev(rowIdx, 0)) == KAT_M::one()) {
492  Kokkos::atomic_add(&numDiagsEqualToOne_dev(), 1);
493  }
494  });
495  }
496  if (useAverageAbsDiagVal) {
497  Teuchos::TimeMonitor MMM = *Teuchos::TimeMonitor::getNewTimer("GetLumpedMatrixDiagonal: useAverageAbsDiagVal");
498  typename Kokkos::View<mag_type, execution_space>::host_mirror_type avgAbsDiagVal = Kokkos::create_mirror_view(avgAbsDiagVal_dev);
499  Kokkos::deep_copy(avgAbsDiagVal, avgAbsDiagVal_dev);
500  int numDiagsEqualToOne;
501  Kokkos::deep_copy(numDiagsEqualToOne, numDiagsEqualToOne_dev);
502 
503  tol = TST::magnitude(100 * Teuchos::ScalarTraits<Scalar>::eps()) * (avgAbsDiagVal() - numDiagsEqualToOne) / (rowMap->getLocalNumElements() - numDiagsEqualToOne);
504  }
505 
506  {
507  Teuchos::TimeMonitor MMM = *Teuchos::TimeMonitor::getNewTimer("ComputeLumpedDiagonalInverse: parallel_for (doReciprocal)");
508  Kokkos::parallel_for(
509  "ComputeLumpedDiagonalInverse", my_policy,
510  KOKKOS_LAMBDA(const int rowIdx) {
511  if (replaceSingleEntryRowWithZero && nnzPerRow(rowIdx) <= 1) {
512  diag_dev(rowIdx, 0) = KAT_S::zero();
513  } else if ((diag_dev(rowIdx, 0) != KAT_S::zero()) && (KAT_S::magnitude(diag_dev(rowIdx, 0)) < KAT_S::magnitude(2 * regSum(rowIdx)))) {
514  diag_dev(rowIdx, 0) = KAT_S::one() / KAT_S::magnitude(2 * regSum(rowIdx));
515  } else {
516  if (KAT_S::magnitude(diag_dev(rowIdx, 0)) > tol) {
517  diag_dev(rowIdx, 0) = KAT_S::one() / diag_dev(rowIdx, 0);
518  } else {
519  diag_dev(rowIdx, 0) = valReplacement_dev;
520  }
521  }
522  });
523  }
524 
525  } else {
526  Teuchos::TimeMonitor MMM = *Teuchos::TimeMonitor::getNewTimer("GetLumpedMatrixDiagonal: parallel_for");
527  Kokkos::parallel_for(
528  "GetLumpedMatrixDiagonal", my_policy,
529  KOKKOS_LAMBDA(const int rowIdx) {
530  diag_dev(rowIdx, 0) = KAT_S::zero();
531  for (size_type entryIdx = local_mat_dev.graph.row_map(rowIdx);
532  entryIdx < local_mat_dev.graph.row_map(rowIdx + 1);
533  ++entryIdx) {
534  diag_dev(rowIdx, 0) += KAT_S::magnitude(local_mat_dev.values(entryIdx));
535  }
536  });
537  }
538  } else {
539  // Implement using Teuchos
540  Teuchos::TimeMonitor MMM = *Teuchos::TimeMonitor::getNewTimer("UtilitiesBase: GetLumpedMatrixDiagonal: (Teuchos implementation)");
541  ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
542  Teuchos::Array<Scalar> regSum(diag->getLocalLength());
543  Teuchos::ArrayView<const LocalOrdinal> cols;
544  Teuchos::ArrayView<const Scalar> vals;
545 
546  std::vector<int> nnzPerRow(rowMap->getLocalNumElements());
547 
548  // FIXME 2021-10-22 JHU If this is called with doReciprocal=false, what should the correct behavior be? Currently,
549  // FIXME 2021-10-22 JHU the diagonal entry is set to be the sum of the absolute values of the row entries.
550 
551  const Magnitude zeroMagn = TST::magnitude(zero);
552  Magnitude avgAbsDiagVal = TST::magnitude(zero);
553  int numDiagsEqualToOne = 0;
554  for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
555  nnzPerRow[i] = 0;
556  rcpA->getLocalRowView(i, cols, vals);
557  diagVals[i] = zero;
558  for (LocalOrdinal j = 0; j < cols.size(); ++j) {
559  regSum[i] += vals[j];
560  const Magnitude rowEntryMagn = TST::magnitude(vals[j]);
561  if (rowEntryMagn > zeroMagn)
562  nnzPerRow[i]++;
563  diagVals[i] += rowEntryMagn;
564  if (static_cast<size_t>(cols[j]) == i)
565  avgAbsDiagVal += rowEntryMagn;
566  }
567  if (nnzPerRow[i] == 1 && TST::magnitude(diagVals[i]) == 1.)
568  numDiagsEqualToOne++;
569  }
570  if (useAverageAbsDiagVal)
571  tol = TST::magnitude(100 * Teuchos::ScalarTraits<Scalar>::eps()) * (avgAbsDiagVal - numDiagsEqualToOne) / (rowMap->getLocalNumElements() - numDiagsEqualToOne);
572  if (doReciprocal) {
573  for (size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
574  if (replaceSingleEntryRowWithZero && nnzPerRow[i] <= static_cast<int>(1))
575  diagVals[i] = zero;
576  else if ((diagVals[i] != zero) && (TST::magnitude(diagVals[i]) < TST::magnitude(two * regSum[i])))
577  diagVals[i] = one / TST::magnitude((two * regSum[i]));
578  else {
579  if (TST::magnitude(diagVals[i]) > tol)
580  diagVals[i] = one / diagVals[i];
581  else {
582  diagVals[i] = valReplacement;
583  }
584  }
585  }
586  }
587  }
588  } else {
589  TEUCHOS_TEST_FOR_EXCEPTION(doReciprocal, Xpetra::Exceptions::RuntimeError,
590  "UtilitiesBase::GetLumpedMatrixDiagonal(): extracting reciprocal of diagonal of a blocked matrix is not supported");
591  diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(bA->getRangeMapExtractor()->getFullMap(), true);
592 
593  for (size_t row = 0; row < bA->Rows(); ++row) {
594  for (size_t col = 0; col < bA->Cols(); ++col) {
595  if (!bA->getMatrix(row, col).is_null()) {
596  // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
597  bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(bA->getMatrix(row, col)) == Teuchos::null);
598  RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag, row, bThyraMode);
599  RCP<const Vector> dd = GetLumpedMatrixDiagonal(*(bA->getMatrix(row, col)));
600  ddtemp->update(Teuchos::as<Scalar>(1.0), *dd, Teuchos::as<Scalar>(1.0));
601  bA->getRangeMapExtractor()->InsertVector(ddtemp, row, diag, bThyraMode);
602  }
603  }
604  }
605  }
606 
607  return diag;
608 }
609 
610 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
611 Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
612 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
613  GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A) {
614  // Get/Create distributed objects
615  RCP<const Map> rowMap = A.getRowMap();
616  auto diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, false);
617 
618  // Implement using Kokkos
619  using local_vector_type = typename Vector::dual_view_type::t_dev_um;
620  using local_matrix_type = typename Matrix::local_matrix_type;
621  using execution_space = typename local_vector_type::execution_space;
622  using values_type = typename local_matrix_type::values_type;
623  using scalar_type = typename values_type::non_const_value_type;
624 #if KOKKOS_VERSION >= 40799
625  using KAT_S = typename KokkosKernels::ArithTraits<scalar_type>;
626 #else
627  using KAT_S = typename Kokkos::ArithTraits<scalar_type>;
628 #endif
629 
630  auto diag_dev = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
631  auto local_mat_dev = A.getLocalMatrixDevice();
632  Kokkos::RangePolicy<execution_space, int> my_policy(0, static_cast<int>(diag_dev.extent(0)));
633 
634  Kokkos::parallel_for(
635  "GetMatrixMaxMinusOffDiagonal", my_policy,
636  KOKKOS_LAMBDA(const LocalOrdinal rowIdx) {
637  auto mymax = KAT_S::zero();
638  auto row = local_mat_dev.rowConst(rowIdx);
639  for (LocalOrdinal entryIdx = 0; entryIdx < row.length; ++entryIdx) {
640  if (rowIdx != row.colidx(entryIdx)) {
641  if (KAT_S::real(mymax) < -KAT_S::real(row.value(entryIdx)))
642  mymax = -KAT_S::real(row.value(entryIdx));
643  }
644  }
645  diag_dev(rowIdx, 0) = mymax;
646  });
647 
648  return diag;
649 }
650 
651 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
652 Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
653 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
654  GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber) {
655  TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()), std::runtime_error, "GetMatrixMaxMinusOffDiagonal: BlockNumber must match's A's column map.");
656 
657  // Get/Create distributed objects
658  RCP<const Map> rowMap = A.getRowMap();
659  auto diag = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(rowMap, false);
660 
661  // Implement using Kokkos
662  using local_vector_type = typename Vector::dual_view_type::t_dev_um;
663  using local_matrix_type = typename Matrix::local_matrix_type;
664  using execution_space = typename local_vector_type::execution_space;
665  using values_type = typename local_matrix_type::values_type;
666  using scalar_type = typename values_type::non_const_value_type;
667 #if KOKKOS_VERSION >= 40799
668  using KAT_S = typename KokkosKernels::ArithTraits<scalar_type>;
669 #else
670  using KAT_S = typename Kokkos::ArithTraits<scalar_type>;
671 #endif
672 
673  auto diag_dev = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
674  auto local_mat_dev = A.getLocalMatrixDevice();
675  auto local_block_dev = BlockNumber.getLocalViewDevice(Xpetra::Access::ReadOnly);
676  Kokkos::RangePolicy<execution_space, int> my_policy(0, static_cast<int>(diag_dev.extent(0)));
677 
678  Kokkos::parallel_for(
679  "GetMatrixMaxMinusOffDiagonal", my_policy,
680  KOKKOS_LAMBDA(const LocalOrdinal rowIdx) {
681  auto mymax = KAT_S::zero();
682  auto row = local_mat_dev.row(rowIdx);
683  for (LocalOrdinal entryIdx = 0; entryIdx < row.length; ++entryIdx) {
684  if ((rowIdx != row.colidx(entryIdx)) && (local_block_dev(rowIdx, 0) == local_block_dev(row.colidx(entryIdx), 0))) {
685  if (KAT_S::real(mymax) < -KAT_S::real(row.value(entryIdx)))
686  mymax = -KAT_S::real(row.value(entryIdx));
687  }
688  }
689  diag_dev(rowIdx, 0) = mymax;
690  });
691 
692  return diag;
693 }
694 
695 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
696 Teuchos::RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
697 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
698  GetInverse(Teuchos::RCP<const Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> v, typename Teuchos::ScalarTraits<Scalar>::magnitudeType tol, Scalar valReplacement) {
699  RCP<Vector> ret = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(v->getMap(), true);
700 
701  // check whether input vector "v" is a BlockedVector
702  RCP<const BlockedVector> bv = Teuchos::rcp_dynamic_cast<const BlockedVector>(v);
703  if (bv.is_null() == false) {
704  RCP<BlockedVector> bret = Teuchos::rcp_dynamic_cast<BlockedVector>(ret);
705  TEUCHOS_TEST_FOR_EXCEPTION(bret.is_null() == true, MueLu::Exceptions::RuntimeError, "MueLu::UtilitiesBase::GetInverse: return vector should be of type BlockedVector");
706  RCP<const BlockedMap> bmap = bv->getBlockedMap();
707  for (size_t r = 0; r < bmap->getNumMaps(); ++r) {
708  RCP<const MultiVector> submvec = bv->getMultiVector(r, bmap->getThyraMode());
709  RCP<const Vector> subvec = submvec->getVector(0);
710  RCP<Vector> subvecinf = MueLu::UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetInverse(subvec, tol, valReplacement);
711  bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
712  }
713  return ret;
714  }
715 
716  // v is an {Epetra,Tpetra}Vector: work with the underlying raw data
717  ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
718  ArrayRCP<const Scalar> inputVals = v->getData(0);
719  for (size_t i = 0; i < v->getMap()->getLocalNumElements(); ++i) {
720  if (Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
721  retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
722  else
723  retVals[i] = valReplacement;
724  }
725  return ret;
726 }
727 
728 // template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
729 // RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
730 // UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
731 // GetMatrixOverlappedDiagonal(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> & A) {
732 // RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
733 
734 // // Undo block map (if we have one)
735 // RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
736 // if(!browMap.is_null()) rowMap = browMap->getMap();
737 
738 // RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
739 // try {
740 // const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
741 // if (crsOp == NULL) {
742 // throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
743 // }
744 // Teuchos::ArrayRCP<size_t> offsets;
745 // crsOp->getLocalDiagOffsets(offsets);
746 // crsOp->getLocalDiagCopy(*localDiag,offsets());
747 // }
748 // catch (...) {
749 // ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
750 // Teuchos::ArrayRCP<Scalar> diagVals = GetMatrixDiagonal(A);
751 // for (LocalOrdinal i = 0; i < localDiagVals.size(); i++)
752 // localDiagVals[i] = diagVals[i];
753 // localDiagVals = diagVals = null;
754 // }
755 
756 // RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
757 // RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
758 // importer = A.getCrsGraph()->getImporter();
759 // if (importer == Teuchos::null) {
760 // importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
761 // }
762 // diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
763 // return diagonal;
764 // }
765 
766 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
767 RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
768 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
769  GetMatrixOverlappedDiagonal(const Matrix& A) {
770  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
771  RCP<Vector> localDiag = GetMatrixDiagonal(A);
772  RCP<const Import> importer = A.getCrsGraph()->getImporter();
773  if (importer.is_null() && !rowMap->isSameAs(*colMap)) {
774  importer = ImportFactory::Build(rowMap, colMap);
775  }
776  if (!importer.is_null()) {
777  RCP<Vector> diagonal = VectorFactory::Build(colMap);
778  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
779  return diagonal;
780  } else
781  return localDiag;
782 }
783 
784 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
785 RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
786 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
787  GetMatrixOverlappedDeletedRowsum(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A) {
788  using STS = typename Teuchos::ScalarTraits<SC>;
789 
790  // Undo block map (if we have one)
791  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
792  RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
793  if (!browMap.is_null()) rowMap = browMap->getMap();
794 
795  RCP<Vector> local = Xpetra::VectorFactory<SC, LO, GO, Node>::Build(rowMap);
796  RCP<Vector> ghosted = Xpetra::VectorFactory<SC, LO, GO, Node>::Build(colMap, true);
797  ArrayRCP<SC> localVals = local->getDataNonConst(0);
798 
799  for (LO row = 0; row < static_cast<LO>(A.getRowMap()->getLocalNumElements()); ++row) {
800  size_t nnz = A.getNumEntriesInLocalRow(row);
801  ArrayView<const LO> indices;
802  ArrayView<const SC> vals;
803  A.getLocalRowView(row, indices, vals);
804 
805  SC si = STS::zero();
806 
807  for (LO colID = 0; colID < static_cast<LO>(nnz); colID++) {
808  if (indices[colID] != row) {
809  si += vals[colID];
810  }
811  }
812  localVals[row] = si;
813  }
814 
815  RCP<const Xpetra::Import<LO, GO, Node>> importer;
816  importer = A.getCrsGraph()->getImporter();
817  if (importer == Teuchos::null) {
818  importer = Xpetra::ImportFactory<LO, GO, Node>::Build(rowMap, colMap);
819  }
820  ghosted->doImport(*local, *(importer), Xpetra::INSERT);
821  return ghosted;
822 }
823 
824 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
825 RCP<Xpetra::Vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node>>
826 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
827  GetMatrixOverlappedAbsDeletedRowsum(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A) {
828  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
829  using STS = typename Teuchos::ScalarTraits<Scalar>;
830  using MTS = typename Teuchos::ScalarTraits<Magnitude>;
831  using MT = Magnitude;
832  using RealValuedVector = Xpetra::Vector<MT, LO, GO, Node>;
833 
834  // Undo block map (if we have one)
835  RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
836  if (!browMap.is_null()) rowMap = browMap->getMap();
837 
838  RCP<RealValuedVector> local = Xpetra::VectorFactory<MT, LO, GO, Node>::Build(rowMap);
839  RCP<RealValuedVector> ghosted = Xpetra::VectorFactory<MT, LO, GO, Node>::Build(colMap, true);
840  ArrayRCP<MT> localVals = local->getDataNonConst(0);
841 
842  for (LO rowIdx = 0; rowIdx < static_cast<LO>(A.getRowMap()->getLocalNumElements()); ++rowIdx) {
843  size_t nnz = A.getNumEntriesInLocalRow(rowIdx);
844  ArrayView<const LO> indices;
845  ArrayView<const SC> vals;
846  A.getLocalRowView(rowIdx, indices, vals);
847 
848  MT si = MTS::zero();
849 
850  for (LO colID = 0; colID < static_cast<LO>(nnz); ++colID) {
851  if (indices[colID] != rowIdx) {
852  si += STS::magnitude(vals[colID]);
853  }
854  }
855  localVals[rowIdx] = si;
856  }
857 
858  RCP<const Xpetra::Import<LO, GO, Node>> importer;
859  importer = A.getCrsGraph()->getImporter();
860  if (importer == Teuchos::null) {
861  importer = Xpetra::ImportFactory<LO, GO, Node>::Build(rowMap, colMap);
862  }
863  ghosted->doImport(*local, *(importer), Xpetra::INSERT);
864  return ghosted;
865 }
866 
867 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
868 GlobalOrdinal
869 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
870  CountNegativeDiagonalEntries(const Matrix& A) {
871  using local_matrix_type = typename Matrix::local_matrix_type;
872  using execution_space = typename local_matrix_type::execution_space;
873 #if KOKKOS_VERSION >= 40799
874  using KAT_S = typename KokkosKernels::ArithTraits<typename local_matrix_type::value_type>;
875 #else
876  using KAT_S = typename Kokkos::ArithTraits<typename local_matrix_type::value_type>;
877 #endif
878 
879  auto local_mat_dev = A.getLocalMatrixDevice();
880  Kokkos::RangePolicy<execution_space, int> my_policy(0, static_cast<int>(local_mat_dev.numRows()));
881  GlobalOrdinal count_l = 0, count_g = 0;
882 
883  Kokkos::parallel_reduce(
884  "CountNegativeDiagonalEntries", my_policy,
885  KOKKOS_LAMBDA(const LocalOrdinal rowIdx, GlobalOrdinal& sum) {
886  auto row = local_mat_dev.row(rowIdx);
887  for (LocalOrdinal entryIdx = 0; entryIdx < row.length; ++entryIdx) {
888  if (rowIdx == row.colidx(entryIdx) && KAT_S::real(row.value(entryIdx)) < 0)
889  sum++;
890  }
891  },
892  count_l);
893 
894  MueLu_sumAll(A.getRowMap()->getComm(), count_l, count_g);
895  return count_g;
896 }
897 
898 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
899 Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>
900 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
901  ResidualNorm(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS) {
902  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
903  const size_t numVecs = X.getNumVectors();
904  RCP<MultiVector> RES = Residual(Op, X, RHS);
905  Teuchos::Array<Magnitude> norms(numVecs);
906  RES->norm2(norms);
907  return norms;
908 }
909 
910 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
911 Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>
912 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
913  ResidualNorm(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector& Resid) {
914  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
915  const size_t numVecs = X.getNumVectors();
916  Residual(Op, X, RHS, Resid);
917  Teuchos::Array<Magnitude> norms(numVecs);
918  Resid.norm2(norms);
919  return norms;
920 }
921 
922 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
923 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
924 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
925  Residual(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS) {
926  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
927  const size_t numVecs = X.getNumVectors();
928  // TODO Op.getRangeMap should return a BlockedMap if it is a BlockedCrsOperator
929  RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(RHS.getMap(), numVecs, false); // no need to initialize to zero
930  Op.residual(X, RHS, *RES);
931  return RES;
932 }
933 
934 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
935 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
936  Residual(const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector& Resid) {
937  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides");
938  TEUCHOS_TEST_FOR_EXCEPTION(Resid.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of residual vectors != number of right-hand sides");
939  Op.residual(X, RHS, Resid);
940 }
941 
942 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
943 Scalar
944 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
945  PowerMethod(const Matrix& A, bool scaleByDiag,
946  LocalOrdinal niters, typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance, bool verbose, unsigned int seed) {
947  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
948  "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
949 
950  // power iteration
951  RCP<Vector> diagInvVec;
952  if (scaleByDiag) {
953  diagInvVec = GetMatrixDiagonalInverse(A);
954  }
955 
956  Scalar lambda = PowerMethod(A, diagInvVec, niters, tolerance, verbose, seed);
957  return lambda;
958 }
959 
960 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
961 Scalar
962 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
963  PowerMethod(const Matrix& A, const RCP<Vector>& diagInvVec,
964  LocalOrdinal niters, typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance, bool verbose, unsigned int seed) {
965  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
966  "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
967 
968  // Create three vectors, fill z with random numbers
969  RCP<Vector> q = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(A.getDomainMap());
970  RCP<Vector> r = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(A.getRangeMap());
971  RCP<Vector> z = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(A.getRangeMap());
972 
973  z->setSeed(seed); // seed random number generator
974  z->randomize(true); // use Xpetra implementation: -> same results for Epetra and Tpetra
975 
976  Teuchos::Array<Magnitude> norms(1);
977 
978  typedef Teuchos::ScalarTraits<Scalar> STS;
979 
980  const Scalar zero = STS::zero(), one = STS::one();
981 
982  Scalar lambda = zero;
983  Magnitude residual = STS::magnitude(zero);
984 
985  // power iteration
986  for (int iter = 0; iter < niters; ++iter) {
987  z->norm2(norms); // Compute 2-norm of z
988  q->update(one / norms[0], *z, zero); // Set q = z / normz
989  A.apply(*q, *z); // Compute z = A*q
990  if (diagInvVec != Teuchos::null)
991  z->elementWiseMultiply(one, *diagInvVec, *z, zero);
992  lambda = q->dot(*z); // Approximate maximum eigenvalue: lamba = dot(q,z)
993 
994  if (iter % 100 == 0 || iter + 1 == niters) {
995  r->update(1.0, *z, -lambda, *q, zero); // Compute A*q - lambda*q
996  r->norm2(norms);
997  residual = STS::magnitude(norms[0] / lambda);
998  if (verbose) {
999  std::cout << "Iter = " << iter
1000  << " Lambda = " << lambda
1001  << " Residual of A*q - lambda*q = " << residual
1002  << std::endl;
1003  }
1004  }
1005  if (residual < tolerance)
1006  break;
1007  }
1008  return lambda;
1009 }
1010 
1011 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1012 RCP<Teuchos::FancyOStream>
1013 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1014  MakeFancy(std::ostream& os) {
1015  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
1016  return fancy;
1017 }
1018 
1019 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1020 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
1021 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1022  Distance2(const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>>& v, LocalOrdinal i0, LocalOrdinal i1) {
1023  const size_t numVectors = v.size();
1024 
1025  Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
1026  for (size_t j = 0; j < numVectors; j++) {
1027  d += (v[j][i0] - v[j][i1]) * (v[j][i0] - v[j][i1]);
1028  }
1029  return Teuchos::ScalarTraits<Scalar>::magnitude(d);
1030 }
1031 
1032 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1033 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
1034 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1035  Distance2(const Teuchos::ArrayView<double>& weight, const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>>& v, LocalOrdinal i0, LocalOrdinal i1) {
1036  const size_t numVectors = v.size();
1037  using MT = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
1038 
1039  Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
1040  for (size_t j = 0; j < numVectors; j++) {
1041  d += Teuchos::as<MT>(weight[j]) * (v[j][i0] - v[j][i1]) * (v[j][i0] - v[j][i1]);
1042  }
1043  return Teuchos::ScalarTraits<Scalar>::magnitude(d);
1044 }
1045 
1046 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1047 Teuchos::ArrayRCP<const bool>
1048 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1049  DetectDirichletRows(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& tol, bool count_twos_as_dirichlet) {
1050  LocalOrdinal numRows = A.getLocalNumRows();
1051  typedef Teuchos::ScalarTraits<Scalar> STS;
1052  ArrayRCP<bool> boundaryNodes(numRows, true);
1053  if (count_twos_as_dirichlet) {
1054  for (LocalOrdinal row = 0; row < numRows; row++) {
1055  ArrayView<const LocalOrdinal> indices;
1056  ArrayView<const Scalar> vals;
1057  A.getLocalRowView(row, indices, vals);
1058  size_t nnz = A.getNumEntriesInLocalRow(row);
1059  if (nnz > 2) {
1060  size_t col;
1061  for (col = 0; col < nnz; col++)
1062  if ((indices[col] != row) && STS::magnitude(vals[col]) > tol) {
1063  if (!boundaryNodes[row])
1064  break;
1065  boundaryNodes[row] = false;
1066  }
1067  if (col == nnz)
1068  boundaryNodes[row] = true;
1069  }
1070  }
1071  } else {
1072  for (LocalOrdinal row = 0; row < numRows; row++) {
1073  ArrayView<const LocalOrdinal> indices;
1074  ArrayView<const Scalar> vals;
1075  A.getLocalRowView(row, indices, vals);
1076  size_t nnz = A.getNumEntriesInLocalRow(row);
1077  if (nnz > 1)
1078  for (size_t col = 0; col < nnz; col++)
1079  if ((indices[col] != row) && STS::magnitude(vals[col]) > tol) {
1080  boundaryNodes[row] = false;
1081  break;
1082  }
1083  }
1084  }
1085  return boundaryNodes;
1086 }
1087 
1088 template <class CrsMatrix>
1089 KOKKOS_FORCEINLINE_FUNCTION bool isDirichletRow(typename CrsMatrix::ordinal_type rowId,
1090  KokkosSparse::SparseRowViewConst<CrsMatrix>& row,
1091 #if KOKKOS_VERSION >= 40799
1092  const typename KokkosKernels::ArithTraits<typename CrsMatrix::value_type>::magnitudeType& tol,
1093 #else
1094  const typename Kokkos::ArithTraits<typename CrsMatrix::value_type>::magnitudeType& tol,
1095 #endif
1096  const bool count_twos_as_dirichlet) {
1097 #if KOKKOS_VERSION >= 40799
1098  using ATS = KokkosKernels::ArithTraits<typename CrsMatrix::value_type>;
1099 #else
1100  using ATS = Kokkos::ArithTraits<typename CrsMatrix::value_type>;
1101 #endif
1102 
1103  auto length = row.length;
1104  bool boundaryNode = true;
1105 
1106  if (count_twos_as_dirichlet) {
1107  if (length > 2) {
1108  decltype(length) colID = 0;
1109  for (; colID < length; colID++)
1110  if ((row.colidx(colID) != rowId) &&
1111  (ATS::magnitude(row.value(colID)) > tol)) {
1112  if (!boundaryNode)
1113  break;
1114  boundaryNode = false;
1115  }
1116  if (colID == length)
1117  boundaryNode = true;
1118  }
1119  } else {
1120  for (decltype(length) colID = 0; colID < length; colID++)
1121  if ((row.colidx(colID) != rowId) &&
1122  (ATS::magnitude(row.value(colID)) > tol)) {
1123  boundaryNode = false;
1124  break;
1125  }
1126  }
1127  return boundaryNode;
1128 }
1129 
1130 template <class SC, class LO, class GO, class NO, class memory_space>
1131 Kokkos::View<bool*, memory_space>
1132 DetectDirichletRows_kokkos(const Xpetra::Matrix<SC, LO, GO, NO>& A,
1133  const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol,
1134  const bool count_twos_as_dirichlet) {
1135 #if KOKKOS_VERSION >= 40799
1136  using impl_scalar_type = typename KokkosKernels::ArithTraits<SC>::val_type;
1137 #else
1138  using impl_scalar_type = typename Kokkos::ArithTraits<SC>::val_type;
1139 #endif
1140 #if KOKKOS_VERSION >= 40799
1141  using ATS = KokkosKernels::ArithTraits<impl_scalar_type>;
1142 #else
1143  using ATS = Kokkos::ArithTraits<impl_scalar_type>;
1144 #endif
1145  using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
1146  using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
1147 
1148  Kokkos::View<bool*, typename NO::device_type::memory_space> boundaryNodes;
1149 
1150  if (helpers::isTpetraBlockCrs(A)) {
1151  const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& Am = toTpetraBlock(A);
1152  auto b_graph = Am.getCrsGraph().getLocalGraphDevice();
1153  auto b_rowptr = Am.getCrsGraph().getLocalRowPtrsDevice();
1154  auto values = Am.getValuesDevice();
1155  LO numBlockRows = Am.getLocalNumRows();
1156  const LO stride = Am.getBlockSize() * Am.getBlockSize();
1157 
1158  boundaryNodes = Kokkos::View<bool*, typename NO::device_type::memory_space>(Kokkos::ViewAllocateWithoutInitializing("boundaryNodes"), numBlockRows);
1159 
1160  if (count_twos_as_dirichlet)
1161  throw Exceptions::RuntimeError("BlockCrs does not support counting twos as Dirichlet");
1162 
1163  Kokkos::parallel_for(
1164  "MueLu:Utils::DetectDirichletRowsBlockCrs", range_type(0, numBlockRows),
1165  KOKKOS_LAMBDA(const LO row) {
1166  auto rowView = b_graph.rowConst(row);
1167  auto length = rowView.length;
1168  LO valstart = b_rowptr[row] * stride;
1169 
1170  boundaryNodes(row) = true;
1171  decltype(length) colID = 0;
1172  for (; colID < length; colID++) {
1173  if (rowView.colidx(colID) != row) {
1174  LO current = valstart + colID * stride;
1175  for (LO k = 0; k < stride; k++) {
1176  if (ATS::magnitude(values[current + k]) > tol) {
1177  boundaryNodes(row) = false;
1178  break;
1179  }
1180  }
1181  }
1182  if (boundaryNodes(row) == false)
1183  break;
1184  }
1185  });
1186  } else {
1187  auto localMatrix = A.getLocalMatrixDevice();
1188  LO numRows = A.getLocalNumRows();
1189  boundaryNodes = Kokkos::View<bool*, typename NO::device_type::memory_space>(Kokkos::ViewAllocateWithoutInitializing("boundaryNodes"), numRows);
1190 
1191  Kokkos::parallel_for(
1192  "MueLu:Utils::DetectDirichletRows", range_type(0, numRows),
1193  KOKKOS_LAMBDA(const LO row) {
1194  auto rowView = localMatrix.rowConst(row);
1195  boundaryNodes(row) = isDirichletRow(row, rowView, tol, count_twos_as_dirichlet);
1196  });
1197  }
1198  if constexpr (std::is_same<memory_space, typename NO::device_type::memory_space>::value)
1199  return boundaryNodes;
1200  else {
1201  Kokkos::View<bool*, memory_space> boundaryNodes2(Kokkos::ViewAllocateWithoutInitializing("boundaryNodes"), boundaryNodes.extent(0));
1202  Kokkos::deep_copy(boundaryNodes2, boundaryNodes);
1203  return boundaryNodes2;
1204  }
1205  // CAG: No idea why this is needed to avoid "warning: missing return statement at end of non-void function"
1206  Kokkos::View<bool*, memory_space> dummy("dummy", 0);
1207  return dummy;
1208 }
1209 
1210 template <class SC, class LO, class GO, class NO>
1211 Kokkos::View<bool*, typename NO::device_type::memory_space>
1212 UtilitiesBase<SC, LO, GO, NO>::
1213  DetectDirichletRows_kokkos(const Xpetra::Matrix<SC, LO, GO, NO>& A,
1214  const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol,
1215  const bool count_twos_as_dirichlet) {
1216  return MueLu::DetectDirichletRows_kokkos<SC, LO, GO, NO, typename NO::device_type::memory_space>(A, tol, count_twos_as_dirichlet);
1217 }
1218 
1219 template <class SC, class LO, class GO, class NO>
1220 Kokkos::View<bool*, typename Kokkos::HostSpace>
1221 UtilitiesBase<SC, LO, GO, NO>::
1222  DetectDirichletRows_kokkos_host(const Xpetra::Matrix<SC, LO, GO, NO>& A,
1223  const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol,
1224  const bool count_twos_as_dirichlet) {
1225  return MueLu::DetectDirichletRows_kokkos<SC, LO, GO, NO, typename Kokkos::HostSpace>(A, tol, count_twos_as_dirichlet);
1226 }
1227 
1228 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1229 Teuchos::ArrayRCP<const bool>
1230 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1231  DetectDirichletRowsExt(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, bool& bHasZeroDiagonal, const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& tol) {
1232  // assume that there is no zero diagonal in matrix
1233  bHasZeroDiagonal = false;
1234 
1235  Teuchos::RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(A.getRowMap());
1236  A.getLocalDiagCopy(*diagVec);
1237  Teuchos::ArrayRCP<const Scalar> diagVecData = diagVec->getData(0);
1238 
1239  LocalOrdinal numRows = A.getLocalNumRows();
1240  typedef Teuchos::ScalarTraits<Scalar> STS;
1241  ArrayRCP<bool> boundaryNodes(numRows, false);
1242  for (LocalOrdinal row = 0; row < numRows; row++) {
1243  ArrayView<const LocalOrdinal> indices;
1244  ArrayView<const Scalar> vals;
1245  A.getLocalRowView(row, indices, vals);
1246  size_t nnz = 0; // collect nonzeros in row (excluding the diagonal)
1247  bool bHasDiag = false;
1248  for (decltype(indices.size()) col = 0; col < indices.size(); col++) {
1249  if (indices[col] != row) {
1250  if (STS::magnitude(vals[col] / STS::magnitude(sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col])))) > tol) {
1251  nnz++;
1252  }
1253  } else
1254  bHasDiag = true; // found a diagonal entry
1255  }
1256  if (bHasDiag == false)
1257  bHasZeroDiagonal = true; // we found at least one row without a diagonal
1258  else if (nnz == 0)
1259  boundaryNodes[row] = true;
1260  }
1261  return boundaryNodes;
1262 }
1263 
1264 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1265 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1266  EnforceInitialCondition(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1267  const Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& RHS,
1268  Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& InitialGuess,
1269  const typename Teuchos::ScalarTraits<SC>::magnitudeType& tol,
1270  const bool count_twos_as_dirichlet) {
1271  using range_type = Kokkos::RangePolicy<LO, typename Node::execution_space>;
1272 
1273  LocalOrdinal numRows = A.getLocalNumRows();
1274  LocalOrdinal numVectors = RHS.getNumVectors();
1275  TEUCHOS_ASSERT_EQUALITY(numVectors, Teuchos::as<LocalOrdinal>(InitialGuess.getNumVectors()));
1276 #ifdef MUELU_DEBUG
1277  TEUCHOS_ASSERT(RHS.getMap()->isCompatible(InitialGuess.getMap()));
1278 #endif
1279 
1280  auto lclRHS = RHS.getLocalViewDevice(Xpetra::Access::ReadOnly);
1281  auto lclInitialGuess = InitialGuess.getLocalViewDevice(Xpetra::Access::ReadWrite);
1282  auto lclA = A.getLocalMatrixDevice();
1283 
1284  Kokkos::parallel_for(
1285  "MueLu:Utils::EnforceInitialCondition", range_type(0, numRows),
1286  KOKKOS_LAMBDA(const LO i) {
1287  auto row = lclA.rowConst(i);
1288  if (isDirichletRow(i, row, tol, count_twos_as_dirichlet)) {
1289  for (LocalOrdinal j = 0; j < numVectors; ++j)
1290  lclInitialGuess(i, j) = lclRHS(i, j);
1291  }
1292  });
1293 }
1294 
1295 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1296 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1297  FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
1298  Teuchos::ArrayRCP<bool> nonzeros) {
1299  TEUCHOS_ASSERT(vals.size() == nonzeros.size());
1300  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
1301  const magnitudeType eps = 2.0 * Teuchos::ScalarTraits<magnitudeType>::eps();
1302  for (size_t i = 0; i < static_cast<size_t>(vals.size()); i++) {
1303  nonzeros[i] = (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > eps);
1304  }
1305 }
1306 
1307 // Find Nonzeros in a device view
1308 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1309 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1310  FindNonZeros(const typename Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::dual_view_type::t_dev_const_um vals,
1311  Kokkos::View<bool*, typename Node::device_type> nonzeros) {
1312 #if KOKKOS_VERSION >= 40799
1313  using ATS = KokkosKernels::ArithTraits<Scalar>;
1314 #else
1315  using ATS = Kokkos::ArithTraits<Scalar>;
1316 #endif
1317 #if KOKKOS_VERSION >= 40799
1318  using impl_ATS = KokkosKernels::ArithTraits<typename ATS::val_type>;
1319 #else
1320  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1321 #endif
1322  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1323  TEUCHOS_ASSERT(vals.extent(0) == nonzeros.extent(0));
1324  const typename ATS::magnitudeType eps = 2.0 * impl_ATS::eps();
1325 
1326  Kokkos::parallel_for(
1327  "MueLu:Maxwell1::FindNonZeros", range_type(0, vals.extent(0)),
1328  KOKKOS_LAMBDA(const size_t i) {
1329  nonzeros(i) = (impl_ATS::magnitude(vals(i, 0)) > eps);
1330  });
1331 }
1332 
1333 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1334 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1335  DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1336  const Teuchos::ArrayRCP<bool>& dirichletRows,
1337  Teuchos::ArrayRCP<bool> dirichletCols,
1338  Teuchos::ArrayRCP<bool> dirichletDomain) {
1339  const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1340  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domMap = A.getDomainMap();
1341  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowMap = A.getRowMap();
1342  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap = A.getColMap();
1343  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == rowMap->getLocalNumElements());
1344  TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == colMap->getLocalNumElements());
1345  TEUCHOS_ASSERT(static_cast<size_t>(dirichletDomain.size()) == domMap->getLocalNumElements());
1346  RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> myColsToZero = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(colMap, 1, /*zeroOut=*/true);
1347  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
1348  for (size_t i = 0; i < (size_t)dirichletRows.size(); i++) {
1349  if (dirichletRows[i]) {
1350  ArrayView<const LocalOrdinal> indices;
1351  ArrayView<const Scalar> values;
1352  A.getLocalRowView(i, indices, values);
1353  for (size_t j = 0; j < static_cast<size_t>(indices.size()); j++)
1354  myColsToZero->replaceLocalValue(indices[j], 0, one);
1355  }
1356  }
1357 
1358  RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> globalColsToZero;
1359  RCP<const Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node>> importer = A.getCrsGraph()->getImporter();
1360  if (!importer.is_null()) {
1361  globalColsToZero = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(domMap, 1, /*zeroOut=*/true);
1362  // export to domain map
1363  globalColsToZero->doExport(*myColsToZero, *importer, Xpetra::ADD);
1364  // import to column map
1365  myColsToZero->doImport(*globalColsToZero, *importer, Xpetra::INSERT);
1366  } else
1367  globalColsToZero = myColsToZero;
1368 
1369  FindNonZeros(globalColsToZero->getData(0), dirichletDomain);
1370  FindNonZeros(myColsToZero->getData(0), dirichletCols);
1371 }
1372 
1373 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1374 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1375  DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1376  const Kokkos::View<bool*, typename Node::device_type>& dirichletRows,
1377  Kokkos::View<bool*, typename Node::device_type> dirichletCols,
1378  Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
1379 #if KOKKOS_VERSION >= 40799
1380  using ATS = KokkosKernels::ArithTraits<Scalar>;
1381 #else
1382  using ATS = Kokkos::ArithTraits<Scalar>;
1383 #endif
1384 #if KOKKOS_VERSION >= 40799
1385  using impl_ATS = KokkosKernels::ArithTraits<typename ATS::val_type>;
1386 #else
1387  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1388 #endif
1389  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1390  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domMap = A.getDomainMap();
1391  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowMap = A.getRowMap();
1392  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap = A.getColMap();
1393  TEUCHOS_ASSERT(dirichletRows.extent(0) == rowMap->getLocalNumElements());
1394  TEUCHOS_ASSERT(dirichletCols.extent(0) == colMap->getLocalNumElements());
1395  TEUCHOS_ASSERT(dirichletDomain.extent(0) == domMap->getLocalNumElements());
1396  RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> myColsToZero = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(colMap, /*zeroOut=*/true);
1397  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
1398  auto myColsToZeroView = myColsToZero->getLocalViewDevice(Xpetra::Access::ReadWrite);
1399  auto localMatrix = A.getLocalMatrixDevice();
1400  Kokkos::parallel_for(
1401  "MueLu:Maxwell1::DetectDirichletCols", range_type(0, rowMap->getLocalNumElements()),
1402  KOKKOS_LAMBDA(const LocalOrdinal row) {
1403  if (dirichletRows(row)) {
1404  auto rowView = localMatrix.row(row);
1405  auto length = rowView.length;
1406 
1407  for (decltype(length) colID = 0; colID < length; colID++)
1408  myColsToZeroView(rowView.colidx(colID), 0) = impl_ATS::one();
1409  }
1410  });
1411 
1412  RCP<Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> globalColsToZero;
1413  RCP<const Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node>> importer = A.getCrsGraph()->getImporter();
1414  if (!importer.is_null()) {
1415  globalColsToZero = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(domMap, /*zeroOut=*/true);
1416  // export to domain map
1417  globalColsToZero->doExport(*myColsToZero, *importer, Xpetra::ADD);
1418  // import to column map
1419  myColsToZero->doImport(*globalColsToZero, *importer, Xpetra::INSERT);
1420  } else
1421  globalColsToZero = myColsToZero;
1422  UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::FindNonZeros(globalColsToZero->getLocalViewDevice(Xpetra::Access::ReadOnly), dirichletDomain);
1423  UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::FindNonZeros(myColsToZero->getLocalViewDevice(Xpetra::Access::ReadOnly), dirichletCols);
1424 }
1425 
1426 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1427 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1428  ApplyRowSumCriterion(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
1429  typedef Teuchos::ScalarTraits<Scalar> STS;
1430  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
1431  typedef Teuchos::ScalarTraits<MT> MTS;
1432  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowmap = A.getRowMap();
1433  for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
1434  size_t nnz = A.getNumEntriesInLocalRow(row);
1435  ArrayView<const LocalOrdinal> indices;
1436  ArrayView<const Scalar> vals;
1437  A.getLocalRowView(row, indices, vals);
1438 
1439  Scalar rowsum = STS::zero();
1440  Scalar diagval = STS::zero();
1441 
1442  for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
1443  LocalOrdinal col = indices[colID];
1444  if (row == col)
1445  diagval = vals[colID];
1446  rowsum += vals[colID];
1447  }
1448  // printf("A(%d,:) row_sum(point) = %6.4e\n",row,rowsum);
1449  if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
1450  // printf("Row %d triggers rowsum\n",(int)row);
1451  dirichletRows[row] = true;
1452  }
1453  }
1454 }
1455 
1456 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1457 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1458  ApplyRowSumCriterion(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber, const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
1459  typedef Teuchos::ScalarTraits<Scalar> STS;
1460  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
1461  typedef Teuchos::ScalarTraits<MT> MTS;
1462  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowmap = A.getRowMap();
1463 
1464  TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()), std::runtime_error, "ApplyRowSumCriterion: BlockNumber must match's A's column map.");
1465 
1466  Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
1467  for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
1468  size_t nnz = A.getNumEntriesInLocalRow(row);
1469  ArrayView<const LocalOrdinal> indices;
1470  ArrayView<const Scalar> vals;
1471  A.getLocalRowView(row, indices, vals);
1472 
1473  Scalar rowsum = STS::zero();
1474  Scalar diagval = STS::zero();
1475  for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
1476  LocalOrdinal col = indices[colID];
1477  if (row == col)
1478  diagval = vals[colID];
1479  if (block_id[row] == block_id[col])
1480  rowsum += vals[colID];
1481  }
1482 
1483  // printf("A(%d,:) row_sum(block) = %6.4e\n",row,rowsum);
1484  if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
1485  // printf("Row %d triggers rowsum\n",(int)row);
1486  dirichletRows[row] = true;
1487  }
1488  }
1489 }
1490 
1491 // Applies rowsum criterion
1492 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class memory_space>
1493 void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1494  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
1495  Kokkos::View<bool*, memory_space>& dirichletRows) {
1496  typedef Teuchos::ScalarTraits<Scalar> STS;
1497  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowmap = A.getRowMap();
1498 
1499  auto dirichletRowsHost = Kokkos::create_mirror_view(dirichletRows);
1500  Kokkos::deep_copy(dirichletRowsHost, dirichletRows);
1501 
1502  for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
1503  size_t nnz = A.getNumEntriesInLocalRow(row);
1504  ArrayView<const LocalOrdinal> indices;
1505  ArrayView<const Scalar> vals;
1506  A.getLocalRowView(row, indices, vals);
1507 
1508  Scalar rowsum = STS::zero();
1509  Scalar diagval = STS::zero();
1510  for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
1511  LocalOrdinal col = indices[colID];
1512  if (row == col)
1513  diagval = vals[colID];
1514  rowsum += vals[colID];
1515  }
1516  if (STS::real(rowsum) > STS::magnitude(diagval) * rowSumTol)
1517  dirichletRowsHost(row) = true;
1518  }
1519 
1520  Kokkos::deep_copy(dirichletRows, dirichletRowsHost);
1521 }
1522 
1523 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1524 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1525  ApplyRowSumCriterion(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1526  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
1527  Kokkos::View<bool*, typename Node::device_type::memory_space>& dirichletRows) {
1528  MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node, typename Node::device_type::memory_space>(A, rowSumTol, dirichletRows);
1529 }
1530 
1531 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1532 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1533  ApplyRowSumCriterionHost(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1534  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
1535  Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows) {
1536  MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node, Kokkos::HostSpace>(A, rowSumTol, dirichletRows);
1537 }
1538 
1539 // Applies rowsum criterion
1540 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class memory_space>
1541 void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1542  const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber,
1543  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
1544  Kokkos::View<bool*, memory_space>& dirichletRows) {
1545  typedef Teuchos::ScalarTraits<Scalar> STS;
1546  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowmap = A.getRowMap();
1547 
1548  TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()), std::runtime_error, "ApplyRowSumCriterion: BlockNumber must match's A's column map.");
1549 
1550  auto dirichletRowsHost = Kokkos::create_mirror_view(dirichletRows);
1551  Kokkos::deep_copy(dirichletRowsHost, dirichletRows);
1552 
1553  Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
1554  for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
1555  size_t nnz = A.getNumEntriesInLocalRow(row);
1556  ArrayView<const LocalOrdinal> indices;
1557  ArrayView<const Scalar> vals;
1558  A.getLocalRowView(row, indices, vals);
1559 
1560  Scalar rowsum = STS::zero();
1561  Scalar diagval = STS::zero();
1562  for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
1563  LocalOrdinal col = indices[colID];
1564  if (row == col)
1565  diagval = vals[colID];
1566  if (block_id[row] == block_id[col])
1567  rowsum += vals[colID];
1568  }
1569  if (STS::real(rowsum) > STS::magnitude(diagval) * rowSumTol)
1570  dirichletRowsHost(row) = true;
1571  }
1572 
1573  Kokkos::deep_copy(dirichletRows, dirichletRowsHost);
1574 }
1575 
1576 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1577 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1578  ApplyRowSumCriterion(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1579  const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber,
1580  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
1581  Kokkos::View<bool*, typename Node::device_type::memory_space>& dirichletRows) {
1582  MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node, typename Node::device_type::memory_space>(A, BlockNumber, rowSumTol, dirichletRows);
1583 }
1584 
1585 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1586 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1587  ApplyRowSumCriterionHost(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1588  const Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>& BlockNumber,
1589  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
1590  Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows) {
1591  MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node, Kokkos::HostSpace>(A, BlockNumber, rowSumTol, dirichletRows);
1592 }
1593 
1594 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1595 Teuchos::ArrayRCP<const bool>
1596 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1597  DetectDirichletCols(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
1598  const Teuchos::ArrayRCP<const bool>& dirichletRows) {
1599  Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1600  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1601  Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domMap = A.getDomainMap();
1602  Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> colMap = A.getColMap();
1603  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> myColsToZero = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(colMap, 1);
1604  myColsToZero->putScalar(zero);
1605  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
1606  for (size_t i = 0; i < (size_t)dirichletRows.size(); i++) {
1607  if (dirichletRows[i]) {
1608  Teuchos::ArrayView<const LocalOrdinal> indices;
1609  Teuchos::ArrayView<const Scalar> values;
1610  A.getLocalRowView(i, indices, values);
1611  for (size_t j = 0; j < static_cast<size_t>(indices.size()); j++)
1612  myColsToZero->replaceLocalValue(indices[j], 0, one);
1613  }
1614  }
1615 
1616  Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>> globalColsToZero = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(domMap, 1);
1617  globalColsToZero->putScalar(zero);
1618  Teuchos::RCP<Xpetra::Export<LocalOrdinal, GlobalOrdinal, Node>> exporter = Xpetra::ExportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(colMap, domMap);
1619  // export to domain map
1620  globalColsToZero->doExport(*myColsToZero, *exporter, Xpetra::ADD);
1621  // import to column map
1622  myColsToZero->doImport(*globalColsToZero, *exporter, Xpetra::INSERT);
1623  Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
1624  Teuchos::ArrayRCP<bool> dirichletCols(colMap->getLocalNumElements(), true);
1625  Magnitude eps = Teuchos::ScalarTraits<Magnitude>::eps();
1626  for (size_t i = 0; i < colMap->getLocalNumElements(); i++) {
1627  dirichletCols[i] = Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i]) > 2.0 * eps;
1628  }
1629  return dirichletCols;
1630 }
1631 
1632 template <class SC, class LO, class GO, class NO>
1633 Kokkos::View<bool*, typename NO::device_type>
1634 UtilitiesBase<SC, LO, GO, NO>::
1635  DetectDirichletCols(const Xpetra::Matrix<SC, LO, GO, NO>& A,
1636  const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows) {
1637 #if KOKKOS_VERSION >= 40799
1638  using ATS = KokkosKernels::ArithTraits<SC>;
1639 #else
1640  using ATS = Kokkos::ArithTraits<SC>;
1641 #endif
1642 #if KOKKOS_VERSION >= 40799
1643  using impl_ATS = KokkosKernels::ArithTraits<typename ATS::val_type>;
1644 #else
1645  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1646 #endif
1647  using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
1648 
1649  SC zero = ATS::zero();
1650 
1651  auto localMatrix = A.getLocalMatrixDevice();
1652  LO numRows = A.getLocalNumRows();
1653 
1654  Teuchos::RCP<const Xpetra::Map<LO, GO, NO>> domMap = A.getDomainMap();
1655  Teuchos::RCP<const Xpetra::Map<LO, GO, NO>> colMap = A.getColMap();
1656  Teuchos::RCP<Xpetra::MultiVector<SC, LO, GO, NO>> myColsToZero = Xpetra::MultiVectorFactory<SC, LO, GO, NO>::Build(colMap, 1);
1657  myColsToZero->putScalar(zero);
1658  auto myColsToZeroView = myColsToZero->getLocalViewDevice(Xpetra::Access::ReadWrite);
1659  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
1660  Kokkos::parallel_for(
1661  "MueLu:Utils::DetectDirichletCols1", range_type(0, numRows),
1662  KOKKOS_LAMBDA(const LO row) {
1663  if (dirichletRows(row)) {
1664  auto rowView = localMatrix.row(row);
1665  auto length = rowView.length;
1666 
1667  for (decltype(length) colID = 0; colID < length; colID++)
1668  myColsToZeroView(rowView.colidx(colID), 0) = impl_ATS::one();
1669  }
1670  });
1671 
1672  Teuchos::RCP<Xpetra::MultiVector<SC, LO, GO, NO>> globalColsToZero = Xpetra::MultiVectorFactory<SC, LO, GO, NO>::Build(domMap, 1);
1673  globalColsToZero->putScalar(zero);
1674  Teuchos::RCP<Xpetra::Export<LO, GO, NO>> exporter = Xpetra::ExportFactory<LO, GO, NO>::Build(colMap, domMap);
1675  // export to domain map
1676  globalColsToZero->doExport(*myColsToZero, *exporter, Xpetra::ADD);
1677  // import to column map
1678  myColsToZero->doImport(*globalColsToZero, *exporter, Xpetra::INSERT);
1679 
1680  auto myCols = myColsToZero->getLocalViewDevice(Xpetra::Access::ReadOnly);
1681  size_t numColEntries = colMap->getLocalNumElements();
1682  Kokkos::View<bool*, typename NO::device_type> dirichletCols(Kokkos::ViewAllocateWithoutInitializing("dirichletCols"), numColEntries);
1683  const typename ATS::magnitudeType eps = 2.0 * ATS::eps();
1684 
1685  Kokkos::parallel_for(
1686  "MueLu:Utils::DetectDirichletCols2", range_type(0, numColEntries),
1687  KOKKOS_LAMBDA(const size_t i) {
1688  dirichletCols(i) = impl_ATS::magnitude(myCols(i, 0)) > eps;
1689  });
1690  return dirichletCols;
1691 }
1692 
1693 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1694 Scalar
1695 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1696  Frobenius(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A, const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& B) {
1697  // We check only row maps. Column may be different. One would hope that they are the same, as we typically
1698  // calculate frobenius norm of the specified sparsity pattern with an updated matrix from the previous step,
1699  // but matrix addition, even when one is submatrix of the other, changes column map (though change may be as
1700  // simple as couple of elements swapped)
1701  TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()), Exceptions::Incompatible, "MueLu::CGSolver::Frobenius: row maps are incompatible");
1702  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(), Exceptions::RuntimeError, "Matrices must be fill completed");
1703 
1704  const Map& AColMap = *A.getColMap();
1705  const Map& BColMap = *B.getColMap();
1706 
1707  Teuchos::ArrayView<const LocalOrdinal> indA, indB;
1708  Teuchos::ArrayView<const Scalar> valA, valB;
1709  size_t nnzA = 0, nnzB = 0;
1710 
1711  // We use a simple algorithm
1712  // for each row we fill valBAll array with the values in the corresponding row of B
1713  // as such, it serves as both sorted array and as storage, so we don't need to do a
1714  // tricky problem: "find a value in the row of B corresponding to the specific GID"
1715  // Once we do that, we translate LID of entries of row of A to LID of B, and multiply
1716  // corresponding entries.
1717  // The algorithm should be reasonably cheap, as it does not sort anything, provided
1718  // that getLocalElement and getGlobalElement functions are reasonably effective. It
1719  // *is* possible that the costs are hidden in those functions, but if maps are close
1720  // to linear maps, we should be fine
1721  Teuchos::Array<Scalar> valBAll(BColMap.getLocalNumElements());
1722 
1723  LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
1724  Scalar zero = Teuchos::ScalarTraits<Scalar>::zero(), f = zero, gf;
1725  size_t numRows = A.getLocalNumRows();
1726  for (size_t i = 0; i < numRows; i++) {
1727  A.getLocalRowView(i, indA, valA);
1728  B.getLocalRowView(i, indB, valB);
1729  nnzA = indA.size();
1730  nnzB = indB.size();
1731 
1732  // Set up array values
1733  for (size_t j = 0; j < nnzB; j++)
1734  valBAll[indB[j]] = valB[j];
1735 
1736  for (size_t j = 0; j < nnzA; j++) {
1737  // The cost of the whole Frobenius dot product function depends on the
1738  // cost of the getLocalElement and getGlobalElement functions here.
1739  LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
1740  if (ind != invalid)
1741  f += valBAll[ind] * valA[j];
1742  }
1743 
1744  // Clean up array values
1745  for (size_t j = 0; j < nnzB; j++)
1746  valBAll[indB[j]] = zero;
1747  }
1748 
1749  MueLu_sumAll(AColMap.getComm(), f, gf);
1750 
1751  return gf;
1752 }
1753 
1754 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1755 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1756  SetRandomSeed(const Teuchos::Comm<int>& comm) {
1757  // Distribute the seeds evenly in [1,maxint-1]. This guarantees nothing
1758  // about where in random number stream we are, but avoids overflow situations
1759  // in parallel when multiplying by a PID. It would be better to use
1760  // a good parallel random number generator.
1761  double one = 1.0;
1762  int maxint = INT_MAX; //= 2^31-1 = 2147483647 for 32-bit integers
1763  int mySeed = Teuchos::as<int>((maxint - 1) * (one - (comm.getRank() + 1) / (comm.getSize() + one)));
1764  if (mySeed < 1 || mySeed == maxint) {
1765  std::ostringstream errStr;
1766  errStr << "Error detected with random seed = " << mySeed << ". It should be in the interval [1,2^31-2].";
1767  throw Exceptions::RuntimeError(errStr.str());
1768  }
1769  std::srand(mySeed);
1770  // For Tpetra, we could use Kokkos' random number generator here.
1771  Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
1772  // Epetra
1773  // MultiVector::Random() -> Epetra_Util::RandomDouble() -> Epetra_Utils::RandomInt()
1774  // Its own random number generator, based on Seed_. Seed_ is initialized in Epetra_Util constructor with std::rand()
1775  // So our setting std::srand() affects that too
1776 }
1777 
1778 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1779 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1780  FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
1781  std::vector<LocalOrdinal>& dirichletRows, bool count_twos_as_dirichlet) {
1782  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
1783  dirichletRows.resize(0);
1784  for (size_t i = 0; i < A->getLocalNumRows(); i++) {
1785  Teuchos::ArrayView<const LocalOrdinal> indices;
1786  Teuchos::ArrayView<const Scalar> values;
1787  A->getLocalRowView(i, indices, values);
1788  int nnz = 0;
1789  for (size_t j = 0; j < (size_t)indices.size(); j++) {
1790  if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > Teuchos::ScalarTraits<MT>::eps()) {
1791  nnz++;
1792  }
1793  }
1794  if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
1795  dirichletRows.push_back(i);
1796  }
1797  }
1798 }
1799 
1800 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1801 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1802  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
1803  const std::vector<LocalOrdinal>& dirichletRows) {
1804  RCP<const Map> Rmap = A->getRowMap();
1805  RCP<const Map> Cmap = A->getColMap();
1806  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1807  Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1808 
1809  for (size_t i = 0; i < dirichletRows.size(); i++) {
1810  GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
1811 
1812  Teuchos::ArrayView<const LocalOrdinal> indices;
1813  Teuchos::ArrayView<const Scalar> values;
1814  A->getLocalRowView(dirichletRows[i], indices, values);
1815  // NOTE: This won't work with fancy node types.
1816  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1817  for (size_t j = 0; j < (size_t)indices.size(); j++) {
1818  if (Cmap->getGlobalElement(indices[j]) == row_gid)
1819  valuesNC[j] = one;
1820  else
1821  valuesNC[j] = zero;
1822  }
1823  }
1824 }
1825 
1826 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1827 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1828  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
1829  const Teuchos::ArrayRCP<const bool>& dirichletRows) {
1830  TEUCHOS_ASSERT(A->isFillComplete());
1831  RCP<const Map> domMap = A->getDomainMap();
1832  RCP<const Map> ranMap = A->getRangeMap();
1833  RCP<const Map> Rmap = A->getRowMap();
1834  RCP<const Map> Cmap = A->getColMap();
1835  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getLocalNumElements());
1836  const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1837  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1838  A->resumeFill();
1839  for (size_t i = 0; i < (size_t)dirichletRows.size(); i++) {
1840  if (dirichletRows[i]) {
1841  GlobalOrdinal row_gid = Rmap->getGlobalElement(i);
1842 
1843  Teuchos::ArrayView<const LocalOrdinal> indices;
1844  Teuchos::ArrayView<const Scalar> values;
1845  A->getLocalRowView(i, indices, values);
1846 
1847  Teuchos::ArrayRCP<Scalar> valuesNC(values.size());
1848  for (size_t j = 0; j < (size_t)indices.size(); j++) {
1849  if (Cmap->getGlobalElement(indices[j]) == row_gid)
1850  valuesNC[j] = one;
1851  else
1852  valuesNC[j] = zero;
1853  }
1854  A->replaceLocalValues(i, indices, valuesNC());
1855  }
1856  }
1857  A->fillComplete(domMap, ranMap);
1858 }
1859 
1860 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1861 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1862  ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
1863  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
1864  TEUCHOS_ASSERT(A->isFillComplete());
1865 #if KOKKOS_VERSION >= 40799
1866  using ATS = KokkosKernels::ArithTraits<Scalar>;
1867 #else
1868  using ATS = Kokkos::ArithTraits<Scalar>;
1869 #endif
1870 #if KOKKOS_VERSION >= 40799
1871  using impl_ATS = KokkosKernels::ArithTraits<typename ATS::val_type>;
1872 #else
1873  using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1874 #endif
1875  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1876 
1877  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> domMap = A->getDomainMap();
1878  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> ranMap = A->getRangeMap();
1879  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> Rmap = A->getRowMap();
1880  RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> Cmap = A->getColMap();
1881 
1882  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getLocalNumElements());
1883 
1884  auto localMatrix = A->getLocalMatrixDevice();
1885  auto localRmap = Rmap->getLocalMap();
1886  auto localCmap = Cmap->getLocalMap();
1887 
1888  Kokkos::parallel_for(
1889  "MueLu::Utils::ApplyOAZ", range_type(0, dirichletRows.extent(0)),
1890  KOKKOS_LAMBDA(const LocalOrdinal row) {
1891  if (dirichletRows(row)) {
1892  auto rowView = localMatrix.row(row);
1893  auto length = rowView.length;
1894  auto row_gid = localRmap.getGlobalElement(row);
1895  auto row_lid = localCmap.getLocalElement(row_gid);
1896 
1897  for (decltype(length) colID = 0; colID < length; colID++)
1898  if (rowView.colidx(colID) == row_lid)
1899  rowView.value(colID) = impl_ATS::one();
1900  else
1901  rowView.value(colID) = impl_ATS::zero();
1902  }
1903  });
1904 }
1905 
1906 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1907 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1908  ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
1909  const std::vector<LocalOrdinal>& dirichletRows,
1910  Scalar replaceWith) {
1911  for (size_t i = 0; i < dirichletRows.size(); i++) {
1912  Teuchos::ArrayView<const LocalOrdinal> indices;
1913  Teuchos::ArrayView<const Scalar> values;
1914  A->getLocalRowView(dirichletRows[i], indices, values);
1915  // NOTE: This won't work with fancy node types.
1916  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1917  for (size_t j = 0; j < (size_t)indices.size(); j++)
1918  valuesNC[j] = replaceWith;
1919  }
1920 }
1921 
1922 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1923 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1924  ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
1925  const Teuchos::ArrayRCP<const bool>& dirichletRows,
1926  Scalar replaceWith) {
1927  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == A->getRowMap()->getLocalNumElements());
1928  for (size_t i = 0; i < (size_t)dirichletRows.size(); i++) {
1929  if (dirichletRows[i]) {
1930  Teuchos::ArrayView<const LocalOrdinal> indices;
1931  Teuchos::ArrayView<const Scalar> values;
1932  A->getLocalRowView(i, indices, values);
1933  // NOTE: This won't work with fancy node types.
1934  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1935  for (size_t j = 0; j < (size_t)indices.size(); j++)
1936  valuesNC[j] = replaceWith;
1937  }
1938  }
1939 }
1940 
1941 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1942 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1943  ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& X,
1944  const Teuchos::ArrayRCP<const bool>& dirichletRows,
1945  Scalar replaceWith) {
1946  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == X->getMap()->getLocalNumElements());
1947  for (size_t i = 0; i < (size_t)dirichletRows.size(); i++) {
1948  if (dirichletRows[i]) {
1949  for (size_t j = 0; j < X->getNumVectors(); j++)
1950  X->replaceLocalValue(i, j, replaceWith);
1951  }
1952  }
1953 }
1954 
1955 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1956 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1957  ZeroDirichletRows(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
1958  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
1959  Scalar replaceWith) {
1960 #if KOKKOS_VERSION >= 40799
1961  using ATS = KokkosKernels::ArithTraits<Scalar>;
1962 #else
1963  using ATS = Kokkos::ArithTraits<Scalar>;
1964 #endif
1965  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1966 
1967  typename ATS::val_type impl_replaceWith = replaceWith;
1968 
1969  auto localMatrix = A->getLocalMatrixDevice();
1970  LocalOrdinal numRows = A->getLocalNumRows();
1971 
1972  Kokkos::parallel_for(
1973  "MueLu:Utils::ZeroDirichletRows", range_type(0, numRows),
1974  KOKKOS_LAMBDA(const LocalOrdinal row) {
1975  if (dirichletRows(row)) {
1976  auto rowView = localMatrix.row(row);
1977  auto length = rowView.length;
1978  for (decltype(length) colID = 0; colID < length; colID++)
1979  rowView.value(colID) = impl_replaceWith;
1980  }
1981  });
1982 }
1983 
1984 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1985 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1986  ZeroDirichletRows(RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& X,
1987  const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
1988  Scalar replaceWith) {
1989 #if KOKKOS_VERSION >= 40799
1990  using ATS = KokkosKernels::ArithTraits<Scalar>;
1991 #else
1992  using ATS = Kokkos::ArithTraits<Scalar>;
1993 #endif
1994  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1995 
1996  typename ATS::val_type impl_replaceWith = replaceWith;
1997 
1998  auto myCols = X->getLocalViewDevice(Xpetra::Access::ReadWrite);
1999  size_t numVecs = X->getNumVectors();
2000  Kokkos::parallel_for(
2001  "MueLu:Utils::ZeroDirichletRows_MV", range_type(0, dirichletRows.size()),
2002  KOKKOS_LAMBDA(const size_t i) {
2003  if (dirichletRows(i)) {
2004  for (size_t j = 0; j < numVecs; j++)
2005  myCols(i, j) = impl_replaceWith;
2006  }
2007  });
2008 }
2009 
2010 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2011 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2012  ZeroDirichletCols(Teuchos::RCP<Matrix>& A,
2013  const Teuchos::ArrayRCP<const bool>& dirichletCols,
2014  Scalar replaceWith) {
2015  TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == A->getColMap()->getLocalNumElements());
2016  for (size_t i = 0; i < A->getLocalNumRows(); i++) {
2017  Teuchos::ArrayView<const LocalOrdinal> indices;
2018  Teuchos::ArrayView<const Scalar> values;
2019  A->getLocalRowView(i, indices, values);
2020  // NOTE: This won't work with fancy node types.
2021  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
2022  for (size_t j = 0; j < static_cast<size_t>(indices.size()); j++)
2023  if (dirichletCols[indices[j]])
2024  valuesNC[j] = replaceWith;
2025  }
2026 }
2027 
2028 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2029 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2030  ZeroDirichletCols(RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
2031  const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols,
2032  Scalar replaceWith) {
2033 #if KOKKOS_VERSION >= 40799
2034  using ATS = KokkosKernels::ArithTraits<Scalar>;
2035 #else
2036  using ATS = Kokkos::ArithTraits<Scalar>;
2037 #endif
2038  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
2039 
2040  typename ATS::val_type impl_replaceWith = replaceWith;
2041 
2042  auto localMatrix = A->getLocalMatrixDevice();
2043  LocalOrdinal numRows = A->getLocalNumRows();
2044 
2045  Kokkos::parallel_for(
2046  "MueLu:Utils::ZeroDirichletCols", range_type(0, numRows),
2047  KOKKOS_LAMBDA(const LocalOrdinal row) {
2048  auto rowView = localMatrix.row(row);
2049  auto length = rowView.length;
2050  for (decltype(length) colID = 0; colID < length; colID++)
2051  if (dirichletCols(rowView.colidx(colID))) {
2052  rowView.value(colID) = impl_replaceWith;
2053  }
2054  });
2055 }
2056 
2057 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2058 void UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2059  FindDirichletRowsAndPropagateToCols(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& A,
2060  Teuchos::RCP<Xpetra::Vector<int, LocalOrdinal, GlobalOrdinal, Node>>& isDirichletRow,
2061  Teuchos::RCP<Xpetra::Vector<int, LocalOrdinal, GlobalOrdinal, Node>>& isDirichletCol) {
2062  // Make sure A's RowMap == DomainMap
2063  if (!A->getRowMap()->isSameAs(*A->getDomainMap())) {
2064  throw std::runtime_error("UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
2065  }
2066  RCP<const Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node>> importer = A->getCrsGraph()->getImporter();
2067  bool has_import = !importer.is_null();
2068 
2069  // Find the Dirichlet rows
2070  std::vector<LocalOrdinal> dirichletRows;
2071  FindDirichletRows(A, dirichletRows);
2072 
2073 #if 0
2074  printf("[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
2075  for(size_t i=0; i<(size_t) dirichletRows.size(); i++)
2076  printf("%d ",dirichletRows[i]);
2077  printf("\n");
2078  fflush(stdout);
2079 #endif
2080  // Allocate all as non-Dirichlet
2081  isDirichletRow = Xpetra::VectorFactory<int, LocalOrdinal, GlobalOrdinal, Node>::Build(A->getRowMap(), true);
2082  isDirichletCol = Xpetra::VectorFactory<int, LocalOrdinal, GlobalOrdinal, Node>::Build(A->getColMap(), true);
2083 
2084  {
2085  Teuchos::ArrayRCP<int> dr_rcp = isDirichletRow->getDataNonConst(0);
2086  Teuchos::ArrayView<int> dr = dr_rcp();
2087  Teuchos::ArrayRCP<int> dc_rcp = isDirichletCol->getDataNonConst(0);
2088  Teuchos::ArrayView<int> dc = dc_rcp();
2089  for (size_t i = 0; i < (size_t)dirichletRows.size(); i++) {
2090  dr[dirichletRows[i]] = 1;
2091  if (!has_import) dc[dirichletRows[i]] = 1;
2092  }
2093  }
2094 
2095  if (has_import)
2096  isDirichletCol->doImport(*isDirichletRow, *importer, Xpetra::CombineMode::ADD);
2097 }
2098 
2099 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2100 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
2101 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2102  ReplaceNonZerosWithOnes(const RCP<Matrix>& original) {
2103 #if KOKKOS_VERSION >= 40799
2104  using ISC = typename KokkosKernels::ArithTraits<Scalar>::val_type;
2105 #else
2106  using ISC = typename Kokkos::ArithTraits<Scalar>::val_type;
2107 #endif
2108  using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
2109  using local_matrix_type = typename CrsMatrix::local_matrix_type;
2110  using values_type = typename local_matrix_type::values_type;
2111 
2112 #if KOKKOS_VERSION >= 40799
2113  const ISC ONE = KokkosKernels::ArithTraits<ISC>::one();
2114 #else
2115  const ISC ONE = Kokkos::ArithTraits<ISC>::one();
2116 #endif
2117 #if KOKKOS_VERSION >= 40799
2118  const ISC ZERO = KokkosKernels::ArithTraits<ISC>::zero();
2119 #else
2120  const ISC ZERO = Kokkos::ArithTraits<ISC>::zero();
2121 #endif
2122 
2123  // Copy the values array of the old matrix to a new array, replacing all the non-zeros with one
2124  auto localMatrix = original->getLocalMatrixDevice();
2125  TEUCHOS_TEST_FOR_EXCEPTION(!original->hasCrsGraph(), Exceptions::RuntimeError, "ReplaceNonZerosWithOnes: Cannot get CrsGraph");
2126  values_type new_values("values", localMatrix.nnz());
2127 
2128  Kokkos::parallel_for(
2129  "ReplaceNonZerosWithOnes", range_type(0, localMatrix.nnz()), KOKKOS_LAMBDA(const size_t i) {
2130  if (localMatrix.values(i) != ZERO)
2131  new_values(i) = ONE;
2132  else
2133  new_values(i) = ZERO;
2134  });
2135 
2136  // Build the new matrix
2137  RCP<Matrix> NewMatrix = Xpetra::MatrixFactory<SC, LO, GO, NO>::Build(original->getCrsGraph(), new_values);
2138  TEUCHOS_TEST_FOR_EXCEPTION(NewMatrix.is_null(), Exceptions::RuntimeError, "ReplaceNonZerosWithOnes: MatrixFactory::Build() did not return matrix");
2139  NewMatrix->fillComplete(original->getDomainMap(), original->getRangeMap());
2140  return NewMatrix;
2141 }
2142 
2143 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2144 RCP<const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>>
2145 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2146  GeneratedBlockedTargetMap(const Xpetra::BlockedMap<LocalOrdinal, GlobalOrdinal, Node>& sourceBlockedMap,
2147  const Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node>& Importer) {
2148  typedef Xpetra::Vector<int, LocalOrdinal, GlobalOrdinal, Node> IntVector;
2149  Xpetra::UnderlyingLib lib = sourceBlockedMap.lib();
2150 
2151  // De-stride the map if we have to (might regret this later)
2152  RCP<const Map> fullMap = sourceBlockedMap.getMap();
2153  RCP<const Map> stridedMap = Teuchos::rcp_dynamic_cast<const Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node>>(fullMap);
2154  if (!stridedMap.is_null()) fullMap = stridedMap->getMap();
2155 
2156  // Initial sanity checking for map compatibil
2157  const size_t numSubMaps = sourceBlockedMap.getNumMaps();
2158  if (!Importer.getSourceMap()->isCompatible(*fullMap))
2159  throw std::runtime_error("GenerateBlockedTargetMap(): Map compatibility error");
2160 
2161  // Build an indicator vector
2162  RCP<IntVector> block_ids = Xpetra::VectorFactory<int, LocalOrdinal, GlobalOrdinal, Node>::Build(fullMap);
2163 
2164  for (size_t i = 0; i < numSubMaps; i++) {
2165  RCP<const Map> map = sourceBlockedMap.getMap(i);
2166 
2167  for (size_t j = 0; j < map->getLocalNumElements(); j++) {
2168  LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
2169  block_ids->replaceLocalValue(jj, (int)i);
2170  }
2171  }
2172 
2173  // Get the block ids for the new map
2174  RCP<const Map> targetMap = Importer.getTargetMap();
2175  RCP<IntVector> new_block_ids = Xpetra::VectorFactory<int, LocalOrdinal, GlobalOrdinal, Node>::Build(targetMap);
2176  new_block_ids->doImport(*block_ids, Importer, Xpetra::CombineMode::ADD);
2177  Teuchos::ArrayRCP<const int> dataRCP = new_block_ids->getData(0);
2178  Teuchos::ArrayView<const int> data = dataRCP();
2179 
2180  // Get the GIDs for each subblock
2181  Teuchos::Array<Teuchos::Array<GlobalOrdinal>> elementsInSubMap(numSubMaps);
2182  for (size_t i = 0; i < targetMap->getLocalNumElements(); i++) {
2183  elementsInSubMap[data[i]].push_back(targetMap->getGlobalElement(i));
2184  }
2185 
2186  // Generate the new submaps
2187  std::vector<RCP<const Map>> subMaps(numSubMaps);
2188  for (size_t i = 0; i < numSubMaps; i++) {
2189  subMaps[i] = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(lib, Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(), elementsInSubMap[i](), targetMap->getIndexBase(), targetMap->getComm());
2190  }
2191 
2192  // Build the BlockedMap
2193  return rcp(new BlockedMap(targetMap, subMaps));
2194 }
2195 
2196 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2197 bool UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2198  MapsAreNested(const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& rowMap, const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& colMap) {
2199  if ((rowMap.lib() == Xpetra::UseTpetra) && (colMap.lib() == Xpetra::UseTpetra)) {
2200  auto tpRowMap = toTpetra(rowMap);
2201  auto tpColMap = toTpetra(colMap);
2202  return tpColMap.isLocallyFitted(tpRowMap);
2203  }
2204 
2205  ArrayView<const GlobalOrdinal> rowElements = rowMap.getLocalElementList();
2206  ArrayView<const GlobalOrdinal> colElements = colMap.getLocalElementList();
2207 
2208  const size_t numElements = rowElements.size();
2209 
2210  if (size_t(colElements.size()) < numElements)
2211  return false;
2212 
2213  bool goodMap = true;
2214  for (size_t i = 0; i < numElements; i++)
2215  if (rowElements[i] != colElements[i]) {
2216  goodMap = false;
2217  break;
2218  }
2219 
2220  return goodMap;
2221 }
2222 
2223 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2224 Teuchos::RCP<Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>
2225 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2226  ReverseCuthillMcKee(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
2227  using local_matrix_type = typename Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_type;
2228  using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
2229  using lno_nnz_view_t = typename local_graph_type::entries_type::non_const_type;
2230  using device = typename local_graph_type::device_type;
2231  using execution_space = typename local_matrix_type::execution_space;
2232  using ordinal_type = typename local_matrix_type::ordinal_type;
2233 
2234  local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
2235 
2236  lno_nnz_view_t rcmOrder = KokkosGraph::Experimental::graph_rcm<device, typename local_graph_type::row_map_type, typename local_graph_type::entries_type, lno_nnz_view_t>(localGraph.row_map, localGraph.entries);
2237 
2238  RCP<Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> retval =
2239  Xpetra::VectorFactory<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>::Build(Op.getRowMap());
2240 
2241  // Copy out and reorder data
2242  auto view1D = Kokkos::subview(retval->getLocalViewDevice(Xpetra::Access::ReadWrite), Kokkos::ALL(), 0);
2243  Kokkos::parallel_for(
2244  "Utilities::ReverseCuthillMcKee",
2245  Kokkos::RangePolicy<ordinal_type, execution_space>(0, localGraph.numRows()),
2246  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
2247  view1D(rcmOrder(rowIdx)) = rowIdx;
2248  });
2249  return retval;
2250 }
2251 
2252 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
2253 Teuchos::RCP<Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>>
2254 UtilitiesBase<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2255  CuthillMcKee(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Op) {
2256  using local_matrix_type = typename Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_type;
2257  using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
2258  using lno_nnz_view_t = typename local_graph_type::entries_type::non_const_type;
2259  using device = typename local_graph_type::device_type;
2260  using execution_space = typename local_matrix_type::execution_space;
2261  using ordinal_type = typename local_matrix_type::ordinal_type;
2262 
2263  local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
2264  LocalOrdinal numRows = localGraph.numRows();
2265 
2266  lno_nnz_view_t rcmOrder = KokkosGraph::Experimental::graph_rcm<device, typename local_graph_type::row_map_type, typename local_graph_type::entries_type, lno_nnz_view_t>(localGraph.row_map, localGraph.entries);
2267 
2268  RCP<Xpetra::Vector<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>> retval =
2269  Xpetra::VectorFactory<LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node>::Build(Op.getRowMap());
2270 
2271  // Copy out data
2272  auto view1D = Kokkos::subview(retval->getLocalViewDevice(Xpetra::Access::ReadWrite), Kokkos::ALL(), 0);
2273  // Since KokkosKernels produced RCM, also reverse the order of the view to get CM
2274  Kokkos::parallel_for(
2275  "Utilities::ReverseCuthillMcKee",
2276  Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
2277  KOKKOS_LAMBDA(const ordinal_type rowIdx) {
2278  view1D(rcmOrder(numRows - 1 - rowIdx)) = rowIdx;
2279  });
2280  return retval;
2281 }
2282 
2283 } // namespace MueLu
2284 
2285 #define MUELU_UTILITIESBASE_SHORT
2286 #endif // MUELU_UTILITIESBASE_DEF_HPP
2287 
2288 // LocalWords: LocalOrdinal
MueLu::UtilitiesBase::FindNonZeros
static void FindNonZeros(const Teuchos::ArrayRCP< const Scalar > vals, Teuchos::ArrayRCP< bool > nonzeros)
Find non-zero values in an ArrayRCP Compares the value to 2 * machine epsilon.
Definition: MueLu_UtilitiesBase_def.hpp:1297
MueLu::isDirichletRow
KOKKOS_FORCEINLINE_FUNCTION bool isDirichletRow(typename CrsMatrix::ordinal_type rowId, KokkosSparse::SparseRowViewConst< CrsMatrix > &row, const typename Kokkos::ArithTraits< typename CrsMatrix::value_type >::magnitudeType &tol, const bool count_twos_as_dirichlet)
Definition: MueLu_UtilitiesBase_def.hpp:1089
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)
MueLu::UtilitiesBase::GetMatrixDiagonal
static RCP< Vector > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
Definition: MueLu_UtilitiesBase_def.hpp:283
MueLu::UtilitiesBase::ApplyOAZToMatrixRows
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const std::vector< LocalOrdinal > &dirichletRows)
Definition: MueLu_UtilitiesBase_def.hpp:1802
MueLu::UtilitiesBase::GetMatrixMaxMinusOffDiagonal
static Teuchos::RCP< Vector > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i=k}(-a_ik), for each for i in the matrix.
Definition: MueLu_UtilitiesBase_def.hpp:613
MueLu::UtilitiesBase::ZeroDirichletCols
static void ZeroDirichletCols(Teuchos::RCP< Matrix > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
Definition: MueLu_UtilitiesBase_def.hpp:2012
MueLu::UtilitiesBase::DetectDirichletColsAndDomains
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomain)
Detects Dirichlet columns &amp; domains from a list of Dirichlet rows.
Definition: MueLu_UtilitiesBase_def.hpp:1335
MueLu_sumAll
#define MueLu_sumAll(rcpComm, in, out)
Definition: MueLu_PerfUtils_decl.hpp:25
Teuchos::Array::append
Array< T > & append(const T &x)
Teuchos::Comm::getSize
virtual int getSize() const =0
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:13
MueLu::ApplyRowSumCriterion
void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, memory_space > &dirichletRows)
Definition: MueLu_UtilitiesBase_def.hpp:1493
Tpetra::BlockCrsMatrix::getBlockSize
virtual LO getBlockSize() const override
MueLu::UtilitiesBase::MapsAreNested
static bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
Definition: MueLu_UtilitiesBase_def.hpp:2198
MueLu::UtilitiesBase::ReverseCuthillMcKee
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Matrix &Op)
Definition: MueLu_UtilitiesBase_def.hpp:2226
MueLu::UtilitiesBase::FindDirichletRowsAndPropagateToCols
static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node >> &isDirichletRow, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node >> &isDirichletCol)
Definition: MueLu_UtilitiesBase_def.hpp:2059
Teuchos::Comm::getRank
virtual int getRank() const =0
MueLu::UtilitiesBase::ApplyRowSumCriterionHost
static void ApplyRowSumCriterionHost(const Matrix &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, Kokkos::HostSpace > &dirichletRows)
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)
Teuchos::ScalarTraits::eps
static magnitudeType eps()
MueLu::removeSmallEntries
Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > removeSmallEntries(Teuchos::RCP< Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType threshold, const bool keepDiagonal)
Definition: MueLu_UtilitiesBase_def.hpp:56
MueLu::UtilitiesBase::PowerMethod
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
MueLu::UtilitiesBase::CountNegativeDiagonalEntries
static GlobalOrdinal CountNegativeDiagonalEntries(const Matrix &A)
Counts the number of negative diagonal entries.
Definition: MueLu_UtilitiesBase_def.hpp:870
MueLu_Exceptions.hpp
MueLu::UtilitiesBase::ApplyRowSumCriterion
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Apply Rowsum Criterion.
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::UtilitiesBase::GetMatrixOverlappedDeletedRowsum
static RCP< Vector > GetMatrixOverlappedDeletedRowsum(const Matrix &A)
Extract Overlapped Matrix Deleted Rowsum.
Definition: MueLu_UtilitiesBase_def.hpp:787
Xpetra::CrsMatrix::getRowMap
virtual const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
Xpetra_CrsGraphFactory.hpp
MueLu::UtilitiesBase::GetThresholdedGraph
static RCP< Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > GetThresholdedGraph(const RCP< Matrix > &A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a graph.
Definition: MueLu_UtilitiesBase_def.hpp:225
Teuchos::ArrayView::size
size_type size() const
Xpetra::Matrix::resumeFill
virtual void resumeFill(const RCP< ParameterList > &params=null)=0
LO
LocalOrdinal LO
MueLu::UtilitiesBase::ZeroDirichletRows
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
Definition: MueLu_UtilitiesBase_def.hpp:1908
MueLu::UtilitiesBase::ReplaceNonZerosWithOnes
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReplaceNonZerosWithOnes(const RCP< Matrix > &original)
Creates a copy of a matrix where the non-zero entries are replaced by ones.
Definition: MueLu_UtilitiesBase_def.hpp:2102
Tpetra::BlockCrsMatrix::getLocalNumRows
size_t getLocalNumRows() const override
MueLu_UtilitiesBase_decl.hpp
Xpetra_CrsMatrixFactory.hpp
Teuchos::ArrayRCP::size
size_type size() const
Xpetra_MatrixMatrix.hpp
MueLu::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: MueLu_Exceptions.hpp:50
MueLu_ConfigDefs.hpp
MueLu::UtilitiesBase::SetRandomSeed
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
Definition: MueLu_UtilitiesBase_def.hpp:1756
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:15
MueLu::UtilitiesBase::CuthillMcKee
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Matrix &Op)
Definition: MueLu_UtilitiesBase_def.hpp:2255
Xpetra_CrsGraph.hpp
Xpetra::Import::getTargetMap
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getTargetMap() const =0
Teuchos::TimeMonitor::getNewTimer
static RCP< Time > getNewTimer(const std::string &name)
Teuchos::ArrayRCP::resize
void resize(const size_type n, const T &val=T())
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
MueLu::UtilitiesBase::Residual
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
Definition: MueLu_UtilitiesBase_def.hpp:925
MueLu::UtilitiesBase::DetectDirichletCols
static Teuchos::ArrayRCP< const bool > DetectDirichletCols(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
Definition: MueLu_UtilitiesBase_def.hpp:1597
Xpetra::Matrix::getLocalRowView
virtual void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const =0
MueLu::UtilitiesBase::FindDirichletRows
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
Definition: MueLu_UtilitiesBase_def.hpp:1780
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:12
MueLu::DetectDirichletRows_kokkos
Kokkos::View< bool *, memory_space > DetectDirichletRows_kokkos(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
Definition: MueLu_UtilitiesBase_def.hpp:1132
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:14
MueLu::UtilitiesBase::GetThresholdedMatrix
static RCP< CrsMatrixWrap > GetThresholdedMatrix(const RCP< Matrix > &Ain, const Magnitude threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a matrix.
Definition: MueLu_UtilitiesBase_def.hpp:169
Xpetra::MapFactory::Build
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
I
#define I(i, j)
MueLu::UtilitiesBase::DetectDirichletRows
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
Definition: MueLu_UtilitiesBase_def.hpp:1049
Xpetra::Matrix::isFillComplete
virtual bool isFillComplete() const =0
MueLu::UtilitiesBase::GetMatrixOverlappedAbsDeletedRowsum
static RCP< Xpetra::Vector< Magnitude, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedAbsDeletedRowsum(const Matrix &A)
Definition: MueLu_UtilitiesBase_def.hpp:827
Xpetra::BlockedMap::getNumMaps
size_t getNumMaps() const
Xpetra::Matrix::replaceLocalValues
virtual void replaceLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)=0
MueLu::UtilitiesBase::GetMatrixDiagonal_arcp
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal_arcp(const Matrix &A)
Extract Matrix Diagonal.
Definition: MueLu_UtilitiesBase_def.hpp:258
Xpetra::MultiVector::getMap
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
MueLu::UtilitiesBase::GeneratedBlockedTargetMap
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
Definition: MueLu_UtilitiesBase_def.hpp:2146
MueLu::UtilitiesBase::GetMatrixOverlappedDiagonal
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
Definition: MueLu_UtilitiesBase_def.hpp:769
Teuchos::ArrayView::getRawPtr
T * getRawPtr() const
MueLu::UtilitiesBase::Magnitude
Teuchos::ScalarTraits< Scalar >::magnitudeType Magnitude
Definition: MueLu_UtilitiesBase_decl.hpp:71
Xpetra::Matrix::getNumEntriesInLocalRow
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
MueLu::UtilitiesBase::GetMatrixDiagonalInverse
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool doLumped=false)
Extract Matrix Diagonal.
Definition: MueLu_UtilitiesBase_def.hpp:319
Teuchos::ScalarTraits::seedrandom
static void seedrandom(unsigned int s)
MueLu::UtilitiesBase::DetectDirichletRowsExt
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows (extended version)
Definition: MueLu_UtilitiesBase_def.hpp:1231
Xpetra::VectorFactory::Build
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
Xpetra::BlockedMap::getMap
const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getMap(size_t i, bool bThyraMode=false) const
Xpetra::BlockedMap::lib
virtual UnderlyingLib lib() const
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
iter
TransListIter iter
toTpetra
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
Teuchos::Array::push_back
void push_back(const value_type &x)
Xpetra::Matrix::getCrsGraph
virtual RCP< const CrsGraph > getCrsGraph() const =0
MueLu::UtilitiesBase::ResidualNorm
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
Definition: MueLu_UtilitiesBase_def.hpp:901
MueLu::UtilitiesBase::EnforceInitialCondition
static void EnforceInitialCondition(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &InitialGuess, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), const bool count_twos_as_dirichlet=false)
Detect Dirichlet rows and copy values from RHS multivector to InitialGuess for Dirichlet rows...
Definition: MueLu_UtilitiesBase_def.hpp:1266
Teuchos::Array::length
int length() const
Teuchos::Array::size
size_type size() const
SC
Scalar SC
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::UtilitiesBase::MakeFancy
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
Definition: MueLu_UtilitiesBase_def.hpp:1014
Xpetra_ExportFactory.hpp
Xpetra::MatrixFactory::Build
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
residual
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:36
MueLu::UtilitiesBase::Frobenius
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)
Frobenius inner product of two matrices.
Definition: MueLu_UtilitiesBase_def.hpp:1696
Xpetra::Import::getSourceMap
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getSourceMap() const =0
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
TEUCHOS_ASSERT_EQUALITY
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
MueLu::UtilitiesBase::GetInverse
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
Definition: MueLu_UtilitiesBase_def.hpp:698
MueLu::UtilitiesBase::DetectDirichletRows_kokkos_host
static Kokkos::View< bool *, typename Kokkos::HostSpace > DetectDirichletRows_kokkos_host(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
Definition: MueLu_UtilitiesBase_def.hpp:1222
Xpetra::MultiVector::getNumVectors
virtual size_t getNumVectors() const =0
MueLu::UtilitiesBase::Distance2
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
Definition: MueLu_UtilitiesBase_def.hpp:1022
MueLu::UtilitiesBase::GetLumpedMatrixDiagonal
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
Definition: MueLu_UtilitiesBase_def.hpp:414
Xpetra_Import.hpp
Xpetra::Map::lib
virtual UnderlyingLib lib() const =0
MueLu::UtilitiesBase::DetectDirichletRows_kokkos
static Kokkos::View< bool *, typename NO::device_type::memory_space > DetectDirichletRows_kokkos(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_UtilitiesBase_def.hpp:1213
MueLu::UtilitiesBase::Crs2Op
static RCP< Matrix > Crs2Op(RCP< CrsMatrix > Op)
Definition: MueLu_UtilitiesBase_def.hpp:45
Xpetra::MultiVector::getData
virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const =0
Xpetra::Operator::getDomainMap
virtual Teuchos::RCP< const Map > getDomainMap() const =0
Xpetra::ImportFactory::Build
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)
Xpetra_ImportFactory.hpp
Teuchos::RCP::is_null
bool is_null() const