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