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