10 #ifndef MUELU_UTILITIESBASE_DEF_HPP
11 #define MUELU_UTILITIESBASE_DEF_HPP
17 #include <Kokkos_Core.hpp>
18 #include <KokkosSparse_CrsMatrix.hpp>
19 #include <KokkosSparse_getDiagCopy.hpp>
21 #include <Xpetra_BlockedVector.hpp>
22 #include <Xpetra_BlockedMap.hpp>
23 #include <Xpetra_BlockedMultiVector.hpp>
31 #include <Xpetra_CrsMatrixWrap.hpp>
32 #include <Xpetra_StridedMap.hpp>
37 #include <KokkosKernels_Handle.hpp>
38 #include <KokkosGraph_RCM.hpp>
42 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
43 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
48 return rcp(
new CrsMatrixWrap(Op));
58 const bool keepDiagonal) {
60 using row_ptr_type =
typename crs_matrix::local_graph_type::row_map_type::non_const_type;
61 using col_idx_type =
typename crs_matrix::local_graph_type::entries_type::non_const_type;
62 using vals_type =
typename crs_matrix::local_matrix_type::values_type;
63 using execution_space =
typename crs_matrix::local_matrix_type::execution_space;
65 using ATS = Kokkos::ArithTraits<Scalar>;
66 using impl_SC =
typename ATS::val_type;
67 using impl_ATS = Kokkos::ArithTraits<impl_SC>;
69 auto lclA = A->getLocalMatrixDevice();
71 auto rowptr = row_ptr_type(
"rowptr", lclA.numRows() + 1);
77 auto lclRowMap = A->
getRowMap()->getLocalMap();
78 auto lclColMap = A->getColMap()->getLocalMap();
79 Kokkos::parallel_scan(
80 "removeSmallEntries::rowptr",
81 Kokkos::RangePolicy<LocalOrdinal, execution_space>(0, lclA.numRows()),
83 auto row = lclA.row(rlid);
84 auto rowInCol = lclColMap.getLocalElement(lclRowMap.getGlobalElement(rlid));
86 if ((impl_ATS::magnitude(row.value(k)) > threshold) || (row.colidx(k) == rowInCol)) {
91 rowptr(rlid + 1) = partial_nnz;
95 idx = col_idx_type(
"idx", nnz);
96 vals = vals_type(
"vals", nnz);
99 "removeSmallEntries::indicesValues",
100 Kokkos::RangePolicy<LocalOrdinal, execution_space>(0, lclA.numRows()),
102 auto row = lclA.row(rlid);
103 auto rowInCol = lclColMap.getLocalElement(lclRowMap.getGlobalElement(rlid));
104 auto I = rowptr(rlid);
106 if ((impl_ATS::magnitude(row.value(k)) > threshold) || (row.colidx(k) == rowInCol)) {
107 idx(
I) = row.colidx(k);
108 vals(
I) = row.value(k);
116 Kokkos::parallel_scan(
117 "removeSmallEntries::rowptr",
118 Kokkos::RangePolicy<LocalOrdinal, execution_space>(0, lclA.numRows()),
120 auto row = lclA.row(rlid);
122 if (impl_ATS::magnitude(row.value(k)) > threshold) {
127 rowptr(rlid + 1) = partial_nnz;
131 idx = col_idx_type(
"idx", nnz);
132 vals = vals_type(
"vals", nnz);
134 Kokkos::parallel_for(
135 "removeSmallEntries::indicesValues",
136 Kokkos::RangePolicy<LocalOrdinal, execution_space>(0, lclA.numRows()),
138 auto row = lclA.row(rlid);
139 auto I = rowptr(rlid);
141 if (impl_ATS::magnitude(row.value(k)) > threshold) {
142 idx(
I) = row.colidx(k);
143 vals(
I) = row.value(k);
152 auto lclNewA =
typename crs_matrix::local_matrix_type(
"thresholdedMatrix", lclA.numRows(), lclA.numCols(), nnz, vals, rowptr, idx);
158 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
159 RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
162 auto crsWrap = rcp_dynamic_cast<CrsMatrixWrap>(Ain);
163 if (!crsWrap.is_null()) {
164 auto crsMat = crsWrap->getCrsMatrix();
166 return rcp_static_cast<CrsMatrixWrap>(filteredMat);
171 RCP<CrsMatrixWrap> Aout =
rcp(
new CrsMatrixWrap(rowmap, expectedNNZperRow <= 0 ? Ain->getGlobalMaxNumRowEntries() : expectedNNZperRow));
173 for (
size_t row = 0; row < Ain->getLocalNumRows(); row++) {
174 size_t nnz = Ain->getNumEntriesInLocalRow(row);
178 Ain->getLocalRowView(row, indices, vals);
184 size_t nNonzeros = 0;
187 LocalOrdinal lclColIdx = colmap->getLocalElement(glbRow);
188 for (
size_t i = 0; i < (size_t)indices.
size(); i++) {
190 indout[nNonzeros] = colmap->getGlobalElement(indices[i]);
191 valout[nNonzeros] = vals[i];
196 for (
size_t i = 0; i < (size_t)indices.
size(); i++) {
198 indout[nNonzeros] = colmap->getGlobalElement(indices[i]);
199 valout[nNonzeros] = vals[i];
205 valout.resize(nNonzeros);
207 Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0, indout.size()), valout.view(0, valout.size()));
209 Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());
214 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
219 RCP<CrsGraph> sparsityPattern = CrsGraphFactory::Build(A->getRowMap(), expectedNNZperRow <= 0 ? A->getGlobalMaxNumRowEntries() : expectedNNZperRow);
221 RCP<Vector> diag = GetMatrixOverlappedDiagonal(*A);
224 for (
size_t row = 0; row < A->getLocalNumRows(); row++) {
227 A->getLocalRowView(row, indices, vals);
229 GlobalOrdinal globalRow = A->getRowMap()->getGlobalElement(row);
230 LocalOrdinal col = A->getColMap()->getLocalElement(globalRow);
232 const Scalar Dk = STS::magnitude(D[col]) > 0.0 ? STS::magnitude(D[col]) : 1.0;
235 for (
size_t i = 0; i < size_t(indices.
size()); i++)
237 if (col == indices[i] || STS::magnitude(STS::squareroot(Dk) * vals[i] * STS::squareroot(Dk)) > STS::magnitude(threshold))
238 indicesNew.
append(A->getColMap()->getGlobalElement(indices[i]));
242 sparsityPattern->fillComplete();
244 return sparsityPattern;
247 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
251 size_t numRows = A.getRowMap()->getLocalNumElements();
255 for (
size_t i = 0; i < numRows; ++i) {
258 for (; j < cols.
size(); ++j) {
259 if (Teuchos::as<size_t>(cols[j]) == i) {
264 if (j == cols.
size()) {
272 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
276 const auto rowMap = A.getRowMap();
279 const CrsMatrixWrap* crsOp =
dynamic_cast<const CrsMatrixWrap*
>(&A);
281 using device_type =
typename CrsGraph::device_type;
282 Kokkos::View<size_t*, device_type> offsets(
"offsets", rowMap->getLocalNumElements());
283 crsOp->getCrsGraph()->getLocalDiagOffsets(offsets);
284 crsOp->getCrsMatrix()->getLocalDiagCopy(*diag, offsets);
286 A.getLocalDiagCopy(*diag);
308 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
314 const bool doLumped) {
327 using local_matrix_type =
typename Matrix::local_matrix_type;
329 using value_type =
typename local_matrix_type::value_type;
330 using values_type =
typename local_matrix_type::values_type;
331 using scalar_type =
typename values_type::non_const_value_type;
332 using ordinal_type =
typename local_matrix_type::ordinal_type;
333 using execution_space =
typename local_matrix_type::execution_space;
340 using KAT = Kokkos::ArithTraits<value_type>;
344 RCP<Vector> diag = VectorFactory::Build(rowMap,
false);
347 local_matrix_type localMatrix = A.getLocalMatrixDevice();
348 auto diagVals = diag->getLocalViewDevice(Xpetra::Access::ReadWrite);
350 ordinal_type numRows = localMatrix.graph.numRows();
352 scalar_type valReplacement_dev = valReplacement;
359 Kokkos::parallel_for(
360 "Utilities::GetMatrixDiagonalInverse",
361 Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
362 KOKKOS_LAMBDA(
const ordinal_type rowIdx) {
363 bool foundDiagEntry =
false;
364 auto myRow = localMatrix.rowConst(rowIdx);
365 for (ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
366 if (myRow.colidx(entryIdx) == rowIdx) {
367 foundDiagEntry =
true;
368 if (KAT::magnitude(myRow.value(entryIdx)) > KAT::magnitude(tol)) {
369 diagVals(rowIdx, 0) = KAT::one() / myRow.value(entryIdx);
371 diagVals(rowIdx, 0) = valReplacement_dev;
377 if (!foundDiagEntry) {
378 diagVals(rowIdx, 0) = KAT::zero();
382 Kokkos::parallel_for(
383 "Utilities::GetMatrixDiagonalInverse",
384 Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
385 KOKKOS_LAMBDA(
const ordinal_type rowIdx) {
386 auto myRow = localMatrix.rowConst(rowIdx);
387 for (ordinal_type entryIdx = 0; entryIdx < myRow.length; ++entryIdx) {
388 diagVals(rowIdx, 0) += KAT::magnitude(myRow.value(entryIdx));
390 if (KAT::magnitude(diagVals(rowIdx, 0)) > KAT::magnitude(tol))
391 diagVals(rowIdx, 0) = KAT::one() / diagVals(rowIdx, 0);
393 diagVals(rowIdx, 0) = valReplacement_dev;
399 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
405 const bool replaceSingleEntryRowWithZero,
406 const bool useAverageAbsDiagVal) {
410 const Scalar zero = TST::zero();
411 const Scalar one = TST::one();
412 const Scalar two = one + one;
418 if (bA == Teuchos::null) {
422 if (rowMap->lib() == Xpetra::UnderlyingLib::UseTpetra) {
425 using local_vector_type =
typename Vector::dual_view_type::t_dev_um;
426 using local_matrix_type =
typename Matrix::local_matrix_type;
427 using execution_space =
typename local_vector_type::execution_space;
430 using values_type =
typename local_matrix_type::values_type;
431 using scalar_type =
typename values_type::non_const_value_type;
432 using mag_type =
typename Kokkos::ArithTraits<scalar_type>::mag_type;
433 using KAT_S =
typename Kokkos::ArithTraits<scalar_type>;
434 using KAT_M =
typename Kokkos::ArithTraits<mag_type>;
435 using size_type =
typename local_matrix_type::non_const_size_type;
437 local_vector_type diag_dev = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
438 local_matrix_type local_mat_dev = rcpA->getLocalMatrixDevice();
439 Kokkos::RangePolicy<execution_space, int> my_policy(0, static_cast<int>(diag_dev.extent(0)));
440 scalar_type valReplacement_dev = valReplacement;
443 Kokkos::View<int*, execution_space> nnzPerRow(
"nnz per rows", diag_dev.extent(0));
444 Kokkos::View<scalar_type*, execution_space> regSum(
"regSum", diag_dev.extent(0));
445 Kokkos::View<mag_type, execution_space> avgAbsDiagVal_dev(
"avgAbsDiagVal");
446 Kokkos::View<int, execution_space> numDiagsEqualToOne_dev(
"numDiagsEqualToOne");
450 Kokkos::parallel_for(
451 "GetLumpedMatrixDiagonal", my_policy,
452 KOKKOS_LAMBDA(
const int rowIdx) {
453 diag_dev(rowIdx, 0) = KAT_S::zero();
454 for (size_type entryIdx = local_mat_dev.graph.row_map(rowIdx);
455 entryIdx < local_mat_dev.graph.row_map(rowIdx + 1);
457 regSum(rowIdx) += local_mat_dev.values(entryIdx);
458 if (KAT_M::zero() < KAT_S::abs(local_mat_dev.values(entryIdx))) {
461 diag_dev(rowIdx, 0) += KAT_S::abs(local_mat_dev.values(entryIdx));
462 if (rowIdx == local_mat_dev.graph.entries(entryIdx)) {
463 Kokkos::atomic_add(&avgAbsDiagVal_dev(), KAT_S::abs(local_mat_dev.values(entryIdx)));
467 if (nnzPerRow(rowIdx) == 1 && KAT_S::magnitude(diag_dev(rowIdx, 0)) == KAT_M::one()) {
468 Kokkos::atomic_add(&numDiagsEqualToOne_dev(), 1);
472 if (useAverageAbsDiagVal) {
474 typename Kokkos::View<mag_type, execution_space>::host_mirror_type avgAbsDiagVal = Kokkos::create_mirror_view(avgAbsDiagVal_dev);
475 Kokkos::deep_copy(avgAbsDiagVal, avgAbsDiagVal_dev);
476 int numDiagsEqualToOne;
477 Kokkos::deep_copy(numDiagsEqualToOne, numDiagsEqualToOne_dev);
484 Kokkos::parallel_for(
485 "ComputeLumpedDiagonalInverse", my_policy,
486 KOKKOS_LAMBDA(
const int rowIdx) {
487 if (replaceSingleEntryRowWithZero && nnzPerRow(rowIdx) <= 1) {
488 diag_dev(rowIdx, 0) = KAT_S::zero();
489 }
else if ((diag_dev(rowIdx, 0) != KAT_S::zero()) && (KAT_S::magnitude(diag_dev(rowIdx, 0)) < KAT_S::magnitude(2 * regSum(rowIdx)))) {
490 diag_dev(rowIdx, 0) = KAT_S::one() / KAT_S::magnitude(2 * regSum(rowIdx));
492 if (KAT_S::magnitude(diag_dev(rowIdx, 0)) > tol) {
493 diag_dev(rowIdx, 0) = KAT_S::one() / diag_dev(rowIdx, 0);
495 diag_dev(rowIdx, 0) = valReplacement_dev;
503 Kokkos::parallel_for(
504 "GetLumpedMatrixDiagonal", my_policy,
505 KOKKOS_LAMBDA(
const int rowIdx) {
506 diag_dev(rowIdx, 0) = KAT_S::zero();
507 for (size_type entryIdx = local_mat_dev.graph.row_map(rowIdx);
508 entryIdx < local_mat_dev.graph.row_map(rowIdx + 1);
510 diag_dev(rowIdx, 0) += KAT_S::magnitude(local_mat_dev.values(entryIdx));
522 std::vector<int> nnzPerRow(rowMap->getLocalNumElements());
527 const Magnitude zeroMagn = TST::magnitude(zero);
528 Magnitude avgAbsDiagVal = TST::magnitude(zero);
529 int numDiagsEqualToOne = 0;
530 for (
size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
532 rcpA->getLocalRowView(i, cols, vals);
535 regSum[i] += vals[j];
536 const Magnitude rowEntryMagn = TST::magnitude(vals[j]);
537 if (rowEntryMagn > zeroMagn)
539 diagVals[i] += rowEntryMagn;
540 if (static_cast<size_t>(cols[j]) == i)
541 avgAbsDiagVal += rowEntryMagn;
543 if (nnzPerRow[i] == 1 && TST::magnitude(diagVals[i]) == 1.)
544 numDiagsEqualToOne++;
546 if (useAverageAbsDiagVal)
549 for (
size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
550 if (replaceSingleEntryRowWithZero && nnzPerRow[i] <= static_cast<int>(1))
552 else if ((diagVals[i] != zero) && (TST::magnitude(diagVals[i]) < TST::magnitude(two * regSum[i])))
553 diagVals[i] = one / TST::magnitude((two * regSum[i]));
555 if (TST::magnitude(diagVals[i]) > tol)
556 diagVals[i] = one / diagVals[i];
558 diagVals[i] = valReplacement;
566 "UtilitiesBase::GetLumpedMatrixDiagonal(): extracting reciprocal of diagonal of a blocked matrix is not supported");
569 for (
size_t row = 0; row < bA->Rows(); ++row) {
570 for (
size_t col = 0; col < bA->Cols(); ++col) {
571 if (!bA->getMatrix(row, col).
is_null()) {
574 RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag, row, bThyraMode);
576 ddtemp->update(Teuchos::as<Scalar>(1.0), *dd, Teuchos::as<Scalar>(1.0));
577 bA->getRangeMapExtractor()->InsertVector(ddtemp, row, diag, bThyraMode);
586 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
595 using local_vector_type =
typename Vector::dual_view_type::t_dev_um;
596 using local_matrix_type =
typename Matrix::local_matrix_type;
597 using execution_space =
typename local_vector_type::execution_space;
598 using values_type =
typename local_matrix_type::values_type;
599 using scalar_type =
typename values_type::non_const_value_type;
600 using KAT_S =
typename Kokkos::ArithTraits<scalar_type>;
602 auto diag_dev = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
603 auto local_mat_dev = A.getLocalMatrixDevice();
604 Kokkos::RangePolicy<execution_space, int> my_policy(0, static_cast<int>(diag_dev.extent(0)));
606 Kokkos::parallel_for(
607 "GetMatrixMaxMinusOffDiagonal", my_policy,
609 auto mymax = KAT_S::zero();
610 auto row = local_mat_dev.rowConst(rowIdx);
611 for (
LocalOrdinal entryIdx = 0; entryIdx < row.length; ++entryIdx) {
612 if (rowIdx != row.colidx(entryIdx)) {
613 if (KAT_S::real(mymax) < -KAT_S::real(row.value(entryIdx)))
614 mymax = -KAT_S::real(row.value(entryIdx));
617 diag_dev(rowIdx, 0) = mymax;
623 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
627 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.
getMap()), std::runtime_error,
"GetMatrixMaxMinusOffDiagonal: BlockNumber must match's A's column map.");
634 using local_vector_type =
typename Vector::dual_view_type::t_dev_um;
635 using local_matrix_type =
typename Matrix::local_matrix_type;
636 using execution_space =
typename local_vector_type::execution_space;
637 using values_type =
typename local_matrix_type::values_type;
638 using scalar_type =
typename values_type::non_const_value_type;
639 using KAT_S =
typename Kokkos::ArithTraits<scalar_type>;
641 auto diag_dev = diag->getLocalViewDevice(Xpetra::Access::OverwriteAll);
642 auto local_mat_dev = A.getLocalMatrixDevice();
643 auto local_block_dev = BlockNumber.getLocalViewDevice(Xpetra::Access::ReadOnly);
644 Kokkos::RangePolicy<execution_space, int> my_policy(0, static_cast<int>(diag_dev.extent(0)));
646 Kokkos::parallel_for(
647 "GetMatrixMaxMinusOffDiagonal", my_policy,
649 auto mymax = KAT_S::zero();
650 auto row = local_mat_dev.row(rowIdx);
651 for (
LocalOrdinal entryIdx = 0; entryIdx < row.length; ++entryIdx) {
652 if ((rowIdx != row.colidx(entryIdx)) && (local_block_dev(rowIdx, 0) == local_block_dev(row.colidx(entryIdx), 0))) {
653 if (KAT_S::real(mymax) < -KAT_S::real(row.value(entryIdx)))
654 mymax = -KAT_S::real(row.value(entryIdx));
657 diag_dev(rowIdx, 0) = mymax;
663 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
671 if (bv.is_null() ==
false) {
675 for (
size_t r = 0; r < bmap->getNumMaps(); ++r) {
679 bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
687 for (
size_t i = 0; i < v->getMap()->getLocalNumElements(); ++i) {
691 retVals[i] = valReplacement;
734 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
741 if (importer.
is_null() && !rowMap->isSameAs(*colMap)) {
742 importer = ImportFactory::Build(rowMap, colMap);
745 RCP<Vector> diagonal = VectorFactory::Build(colMap);
752 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
761 if (!browMap.
is_null()) rowMap = browMap->getMap();
767 for (
LO row = 0; row < static_cast<LO>(A.getRowMap()->getLocalNumElements()); ++row) {
775 for (
LO colID = 0; colID < static_cast<LO>(nnz); colID++) {
776 if (indices[colID] != row) {
785 if (importer == Teuchos::null) {
792 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
804 if (!browMap.
is_null()) rowMap = browMap->getMap();
810 for (
LO rowIdx = 0; rowIdx < static_cast<LO>(A.getRowMap()->getLocalNumElements()); ++rowIdx) {
818 for (
LO colID = 0; colID < static_cast<LO>(nnz); ++colID) {
819 if (indices[colID] != rowIdx) {
820 si += STS::magnitude(vals[colID]);
823 localVals[rowIdx] = si;
828 if (importer == Teuchos::null) {
835 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
840 const size_t numVecs = X.getNumVectors();
847 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
852 const size_t numVecs = X.getNumVectors();
853 Residual(Op, X, RHS, Resid);
859 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
864 const size_t numVecs = X.getNumVectors();
867 Op.residual(X, RHS, *RES);
871 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
876 Op.residual(X, RHS, Resid);
879 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
885 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
890 diagInvVec = GetMatrixDiagonalInverse(A);
893 Scalar lambda = PowerMethod(A, diagInvVec, niters, tolerance, verbose, seed);
897 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
900 PowerMethod(
const Matrix& A,
const RCP<Vector>& diagInvVec,
903 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
917 const Scalar zero = STS::zero(), one = STS::one();
920 Magnitude
residual = STS::magnitude(zero);
923 for (
int iter = 0;
iter < niters; ++
iter) {
925 q->update(one / norms[0], *z, zero);
927 if (diagInvVec != Teuchos::null)
928 z->elementWiseMultiply(one, *diagInvVec, *z, zero);
931 if (iter % 100 == 0 || iter + 1 == niters) {
932 r->update(1.0, *z, -lambda, *q, zero);
934 residual = STS::magnitude(norms[0] / lambda);
936 std::cout <<
"Iter = " <<
iter
937 <<
" Lambda = " << lambda
938 <<
" Residual of A*q - lambda*q = " << residual
942 if (residual < tolerance)
948 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
949 RCP<Teuchos::FancyOStream>
956 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
960 const size_t numVectors = v.size();
963 for (
size_t j = 0; j < numVectors; j++) {
964 d += (v[j][i0] - v[j][i1]) * (v[j][i0] - v[j][i1]);
969 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
973 const size_t numVectors = v.size();
977 for (
size_t j = 0; j < numVectors; j++) {
978 d += Teuchos::as<MT>(weight[j]) * (v[j][i0] - v[j][i1]) * (v[j][i0] - v[j][i1]);
983 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
990 if (count_twos_as_dirichlet) {
998 for (col = 0; col < nnz; col++)
999 if ((indices[col] != row) && STS::magnitude(vals[col]) > tol) {
1000 if (!boundaryNodes[row])
1002 boundaryNodes[row] =
false;
1005 boundaryNodes[row] =
true;
1015 for (
size_t col = 0; col < nnz; col++)
1016 if ((indices[col] != row) && STS::magnitude(vals[col]) > tol) {
1017 boundaryNodes[row] =
false;
1022 return boundaryNodes;
1025 template <
class SC,
class LO,
class GO,
class NO,
class memory_space>
1026 Kokkos::View<bool*, memory_space>
1029 const bool count_twos_as_dirichlet) {
1030 using impl_scalar_type =
typename Kokkos::ArithTraits<SC>::val_type;
1031 using ATS = Kokkos::ArithTraits<impl_scalar_type>;
1032 using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
1035 Kokkos::View<bool*, typename NO::device_type::memory_space> boundaryNodes;
1037 if (helpers::isTpetraBlockCrs(A)) {
1039 auto b_graph = Am.getCrsGraph().getLocalGraphDevice();
1040 auto b_rowptr = Am.getCrsGraph().getLocalRowPtrsDevice();
1045 boundaryNodes = Kokkos::View<bool*, typename NO::device_type::memory_space>(Kokkos::ViewAllocateWithoutInitializing(
"boundaryNodes"), numBlockRows);
1047 if (count_twos_as_dirichlet)
1050 Kokkos::parallel_for(
1051 "MueLu:Utils::DetectDirichletRowsBlockCrs", range_type(0, numBlockRows),
1052 KOKKOS_LAMBDA(
const LO row) {
1053 auto rowView = b_graph.rowConst(row);
1054 auto length = rowView.length;
1055 LO valstart = b_rowptr[row] * stride;
1057 boundaryNodes(row) =
true;
1058 decltype(length) colID = 0;
1059 for (; colID < length; colID++) {
1060 if (rowView.colidx(colID) != row) {
1061 LO current = valstart + colID * stride;
1062 for (
LO k = 0; k < stride; k++) {
1063 if (ATS::magnitude(values[current + k]) > tol) {
1064 boundaryNodes(row) =
false;
1069 if (boundaryNodes(row) ==
false)
1074 auto localMatrix = A.getLocalMatrixDevice();
1075 LO numRows = A.getLocalNumRows();
1076 boundaryNodes = Kokkos::View<bool*, typename NO::device_type::memory_space>(Kokkos::ViewAllocateWithoutInitializing(
"boundaryNodes"), numRows);
1078 if (count_twos_as_dirichlet)
1079 Kokkos::parallel_for(
1080 "MueLu:Utils::DetectDirichletRows_Twos_As_Dirichlet", range_type(0, numRows),
1081 KOKKOS_LAMBDA(
const LO row) {
1082 auto rowView = localMatrix.row(row);
1083 auto length = rowView.length;
1085 boundaryNodes(row) =
true;
1087 decltype(length) colID = 0;
1088 for (; colID < length; colID++)
1089 if ((rowView.colidx(colID) != row) &&
1090 (ATS::magnitude(rowView.value(colID)) > tol)) {
1091 if (!boundaryNodes(row))
1093 boundaryNodes(row) =
false;
1095 if (colID == length)
1096 boundaryNodes(row) =
true;
1100 Kokkos::parallel_for(
1101 "MueLu:Utils::DetectDirichletRows", range_type(0, numRows),
1102 KOKKOS_LAMBDA(
const LO row) {
1103 auto rowView = localMatrix.row(row);
1104 auto length = rowView.length;
1106 boundaryNodes(row) =
true;
1107 for (decltype(length) colID = 0; colID < length; colID++)
1108 if ((rowView.colidx(colID) != row) &&
1109 (ATS::magnitude(rowView.value(colID)) > tol)) {
1110 boundaryNodes(row) =
false;
1115 if constexpr (std::is_same<memory_space, typename NO::device_type::memory_space>::value)
1116 return boundaryNodes;
1118 Kokkos::View<bool*, memory_space> boundaryNodes2(Kokkos::ViewAllocateWithoutInitializing(
"boundaryNodes"), boundaryNodes.extent(0));
1119 Kokkos::deep_copy(boundaryNodes2, boundaryNodes);
1120 return boundaryNodes2;
1123 Kokkos::View<bool*, memory_space> dummy(
"dummy", 0);
1127 template <
class SC,
class LO,
class GO,
class NO>
1128 Kokkos::View<bool*, typename NO::device_type::memory_space>
1132 const bool count_twos_as_dirichlet) {
1133 return MueLu::DetectDirichletRows_kokkos<SC, LO, GO, NO, typename NO::device_type::memory_space>(A, tol, count_twos_as_dirichlet);
1136 template <
class SC,
class LO,
class GO,
class NO>
1137 Kokkos::View<bool*, typename Kokkos::HostSpace>
1141 const bool count_twos_as_dirichlet) {
1142 return MueLu::DetectDirichletRows_kokkos<SC, LO, GO, NO, typename Kokkos::HostSpace>(A, tol, count_twos_as_dirichlet);
1145 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1150 bHasZeroDiagonal =
false;
1164 bool bHasDiag =
false;
1165 for (decltype(indices.
size()) col = 0; col < indices.
size(); col++) {
1166 if (indices[col] != row) {
1167 if (STS::magnitude(vals[col] / STS::magnitude(sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col])))) > tol) {
1173 if (bHasDiag ==
false)
1174 bHasZeroDiagonal =
true;
1176 boundaryNodes[row] =
true;
1178 return boundaryNodes;
1181 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1187 const bool count_twos_as_dirichlet) {
1188 using range_type = Kokkos::RangePolicy<LO, typename Node::execution_space>;
1199 auto lclRHS = RHS.getLocalViewDevice(Xpetra::Access::ReadOnly);
1200 auto lclInitialGuess = InitialGuess.getLocalViewDevice(Xpetra::Access::ReadWrite);
1202 Kokkos::parallel_for(
1203 "MueLu:Utils::EnforceInitialCondition", range_type(0, numRows),
1204 KOKKOS_LAMBDA(
const LO row) {
1205 if (dirichletRows(row)) {
1207 lclInitialGuess(row, j) = lclRHS(row, j);
1212 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1219 for (
size_t i = 0; i < static_cast<size_t>(vals.
size()); i++) {
1225 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1228 Kokkos::View<bool*, typename Node::device_type> nonzeros) {
1229 using ATS = Kokkos::ArithTraits<Scalar>;
1230 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1231 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1233 const typename ATS::magnitudeType eps = 2.0 * impl_ATS::eps();
1235 Kokkos::parallel_for(
1236 "MueLu:Maxwell1::FindNonZeros", range_type(0, vals.extent(0)),
1237 KOKKOS_LAMBDA(
const size_t i) {
1238 nonzeros(i) = (impl_ATS::magnitude(vals(i, 0)) > eps);
1242 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1252 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.
size()) == rowMap->getLocalNumElements());
1253 TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.
size()) == colMap->getLocalNumElements());
1254 TEUCHOS_ASSERT(static_cast<size_t>(dirichletDomain.
size()) == domMap->getLocalNumElements());
1257 for (
size_t i = 0; i < (size_t)dirichletRows.
size(); i++) {
1258 if (dirichletRows[i]) {
1262 for (
size_t j = 0; j < static_cast<size_t>(indices.
size()); j++)
1263 myColsToZero->replaceLocalValue(indices[j], 0, one);
1272 globalColsToZero->doExport(*myColsToZero, *importer,
Xpetra::ADD);
1274 myColsToZero->doImport(*globalColsToZero, *importer,
Xpetra::INSERT);
1276 globalColsToZero = myColsToZero;
1278 FindNonZeros(globalColsToZero->getData(0), dirichletDomain);
1279 FindNonZeros(myColsToZero->getData(0), dirichletCols);
1282 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1285 const Kokkos::View<bool*, typename Node::device_type>& dirichletRows,
1286 Kokkos::View<bool*, typename Node::device_type> dirichletCols,
1287 Kokkos::View<bool*, typename Node::device_type> dirichletDomain) {
1288 using ATS = Kokkos::ArithTraits<Scalar>;
1289 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1290 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1294 TEUCHOS_ASSERT(dirichletRows.extent(0) == rowMap->getLocalNumElements());
1295 TEUCHOS_ASSERT(dirichletCols.extent(0) == colMap->getLocalNumElements());
1296 TEUCHOS_ASSERT(dirichletDomain.extent(0) == domMap->getLocalNumElements());
1299 auto myColsToZeroView = myColsToZero->getLocalViewDevice(Xpetra::Access::ReadWrite);
1300 auto localMatrix = A.getLocalMatrixDevice();
1301 Kokkos::parallel_for(
1302 "MueLu:Maxwell1::DetectDirichletCols", range_type(0, rowMap->getLocalNumElements()),
1304 if (dirichletRows(row)) {
1305 auto rowView = localMatrix.row(row);
1306 auto length = rowView.length;
1308 for (decltype(length) colID = 0; colID < length; colID++)
1309 myColsToZeroView(rowView.colidx(colID), 0) = impl_ATS::one();
1318 globalColsToZero->doExport(*myColsToZero, *importer,
Xpetra::ADD);
1320 myColsToZero->doImport(*globalColsToZero, *importer,
Xpetra::INSERT);
1322 globalColsToZero = myColsToZero;
1327 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1334 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
1340 Scalar rowsum = STS::zero();
1341 Scalar diagval = STS::zero();
1343 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
1346 diagval = vals[colID];
1347 rowsum += vals[colID];
1350 if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
1352 dirichletRows[row] =
true;
1357 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1363 RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> rowmap = A.getRowMap();
1365 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.
getMap()), std::runtime_error,
"ApplyRowSumCriterion: BlockNumber must match's A's column map.");
1368 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
1370 ArrayView<const LocalOrdinal> indices;
1371 ArrayView<const Scalar> vals;
1374 Scalar rowsum = STS::zero();
1375 Scalar diagval = STS::zero();
1376 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
1379 diagval = vals[colID];
1380 if (block_id[row] == block_id[col])
1381 rowsum += vals[colID];
1385 if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
1387 dirichletRows[row] =
true;
1393 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class memory_space>
1396 Kokkos::View<bool*, memory_space>& dirichletRows) {
1400 auto dirichletRowsHost = Kokkos::create_mirror_view(dirichletRows);
1401 Kokkos::deep_copy(dirichletRowsHost, dirichletRows);
1403 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
1409 Scalar rowsum = STS::zero();
1410 Scalar diagval = STS::zero();
1411 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
1414 diagval = vals[colID];
1415 rowsum += vals[colID];
1417 if (STS::real(rowsum) > STS::magnitude(diagval) * rowSumTol)
1418 dirichletRowsHost(row) =
true;
1421 Kokkos::deep_copy(dirichletRows, dirichletRowsHost);
1424 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1428 Kokkos::View<bool*, typename Node::device_type::memory_space>& dirichletRows) {
1429 MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node, typename Node::device_type::memory_space>(A, rowSumTol, dirichletRows);
1432 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1436 Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows) {
1437 MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node, Kokkos::HostSpace>(A, rowSumTol, dirichletRows);
1441 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class memory_space>
1445 Kokkos::View<bool*, memory_space>& dirichletRows) {
1449 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.
getMap()), std::runtime_error,
"ApplyRowSumCriterion: BlockNumber must match's A's column map.");
1451 auto dirichletRowsHost = Kokkos::create_mirror_view(dirichletRows);
1452 Kokkos::deep_copy(dirichletRowsHost, dirichletRows);
1455 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
1461 Scalar rowsum = STS::zero();
1462 Scalar diagval = STS::zero();
1463 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
1466 diagval = vals[colID];
1467 if (block_id[row] == block_id[col])
1468 rowsum += vals[colID];
1470 if (STS::real(rowsum) > STS::magnitude(diagval) * rowSumTol)
1471 dirichletRowsHost(row) =
true;
1474 Kokkos::deep_copy(dirichletRows, dirichletRowsHost);
1477 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1482 Kokkos::View<bool*, typename Node::device_type::memory_space>& dirichletRows) {
1483 MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node, typename Node::device_type::memory_space>(A, BlockNumber, rowSumTol, dirichletRows);
1486 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1491 Kokkos::View<bool*, Kokkos::HostSpace>& dirichletRows) {
1492 MueLu::ApplyRowSumCriterion<Scalar, LocalOrdinal, GlobalOrdinal, Node, Kokkos::HostSpace>(A, BlockNumber, rowSumTol, dirichletRows);
1495 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1505 myColsToZero->putScalar(zero);
1507 for (
size_t i = 0; i < (size_t)dirichletRows.
size(); i++) {
1508 if (dirichletRows[i]) {
1512 for (
size_t j = 0; j < static_cast<size_t>(indices.
size()); j++)
1513 myColsToZero->replaceLocalValue(indices[j], 0, one);
1518 globalColsToZero->putScalar(zero);
1521 globalColsToZero->doExport(*myColsToZero, *exporter,
Xpetra::ADD);
1523 myColsToZero->doImport(*globalColsToZero, *exporter,
Xpetra::INSERT);
1527 for (
size_t i = 0; i < colMap->getLocalNumElements(); i++) {
1530 return dirichletCols;
1533 template <
class SC,
class LO,
class GO,
class NO>
1534 Kokkos::View<bool*, typename NO::device_type>
1537 const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows) {
1538 using ATS = Kokkos::ArithTraits<SC>;
1539 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1540 using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
1542 SC zero = ATS::zero();
1544 auto localMatrix = A.getLocalMatrixDevice();
1545 LO numRows = A.getLocalNumRows();
1550 myColsToZero->putScalar(zero);
1551 auto myColsToZeroView = myColsToZero->getLocalViewDevice(Xpetra::Access::ReadWrite);
1553 Kokkos::parallel_for(
1554 "MueLu:Utils::DetectDirichletCols1", range_type(0, numRows),
1555 KOKKOS_LAMBDA(
const LO row) {
1556 if (dirichletRows(row)) {
1557 auto rowView = localMatrix.row(row);
1558 auto length = rowView.length;
1560 for (decltype(length) colID = 0; colID < length; colID++)
1561 myColsToZeroView(rowView.colidx(colID), 0) = impl_ATS::one();
1566 globalColsToZero->putScalar(zero);
1569 globalColsToZero->doExport(*myColsToZero, *exporter,
Xpetra::ADD);
1571 myColsToZero->doImport(*globalColsToZero, *exporter,
Xpetra::INSERT);
1573 auto myCols = myColsToZero->getLocalViewDevice(Xpetra::Access::ReadOnly);
1574 size_t numColEntries = colMap->getLocalNumElements();
1575 Kokkos::View<bool*, typename NO::device_type> dirichletCols(Kokkos::ViewAllocateWithoutInitializing(
"dirichletCols"), numColEntries);
1576 const typename ATS::magnitudeType eps = 2.0 * ATS::eps();
1578 Kokkos::parallel_for(
1579 "MueLu:Utils::DetectDirichletCols2", range_type(0, numColEntries),
1580 KOKKOS_LAMBDA(
const size_t i) {
1581 dirichletCols(i) = impl_ATS::magnitude(myCols(i, 0)) > eps;
1583 return dirichletCols;
1586 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1597 const Map& AColMap = *A.getColMap();
1598 const Map& BColMap = *B.getColMap();
1602 size_t nnzA = 0, nnzB = 0;
1618 size_t numRows = A.getLocalNumRows();
1619 for (
size_t i = 0; i < numRows; i++) {
1626 for (
size_t j = 0; j < nnzB; j++)
1627 valBAll[indB[j]] = valB[j];
1629 for (
size_t j = 0; j < nnzA; j++) {
1632 LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
1634 f += valBAll[ind] * valA[j];
1638 for (
size_t j = 0; j < nnzB; j++)
1639 valBAll[indB[j]] = zero;
1647 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1655 int maxint = INT_MAX;
1656 int mySeed = Teuchos::as<int>((maxint - 1) * (one - (comm.
getRank() + 1) / (comm.
getSize() + one)));
1657 if (mySeed < 1 || mySeed == maxint) {
1658 std::ostringstream errStr;
1659 errStr <<
"Error detected with random seed = " << mySeed <<
". It should be in the interval [1,2^31-2].";
1671 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1674 std::vector<LocalOrdinal>& dirichletRows,
bool count_twos_as_dirichlet) {
1676 dirichletRows.resize(0);
1677 for (
size_t i = 0; i < A->getLocalNumRows(); i++) {
1682 for (
size_t j = 0; j < (size_t)indices.
size(); j++) {
1687 if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
1688 dirichletRows.push_back(i);
1693 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1696 const std::vector<LocalOrdinal>& dirichletRows) {
1702 for (
size_t i = 0; i < dirichletRows.size(); i++) {
1703 GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
1710 for (
size_t j = 0; j < (size_t)indices.
size(); j++) {
1711 if (Cmap->getGlobalElement(indices[j]) == row_gid)
1719 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1728 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.
size()) == Rmap->getLocalNumElements());
1732 for (
size_t i = 0; i < (size_t)dirichletRows.
size(); i++) {
1733 if (dirichletRows[i]) {
1741 for (
size_t j = 0; j < (size_t)indices.
size(); j++) {
1742 if (Cmap->getGlobalElement(indices[j]) == row_gid)
1753 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1756 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows) {
1758 using ATS = Kokkos::ArithTraits<Scalar>;
1759 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1760 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1767 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getLocalNumElements());
1769 auto localMatrix = A->getLocalMatrixDevice();
1770 auto localRmap = Rmap->getLocalMap();
1771 auto localCmap = Cmap->getLocalMap();
1773 Kokkos::parallel_for(
1774 "MueLu::Utils::ApplyOAZ", range_type(0, dirichletRows.extent(0)),
1776 if (dirichletRows(row)) {
1777 auto rowView = localMatrix.row(row);
1778 auto length = rowView.length;
1779 auto row_gid = localRmap.getGlobalElement(row);
1780 auto row_lid = localCmap.getLocalElement(row_gid);
1782 for (decltype(length) colID = 0; colID < length; colID++)
1783 if (rowView.colidx(colID) == row_lid)
1784 rowView.value(colID) = impl_ATS::one();
1786 rowView.value(colID) = impl_ATS::zero();
1791 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1794 const std::vector<LocalOrdinal>& dirichletRows,
1796 for (
size_t i = 0; i < dirichletRows.size(); i++) {
1802 for (
size_t j = 0; j < (size_t)indices.
size(); j++)
1803 valuesNC[j] = replaceWith;
1807 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1812 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.
size()) == A->getRowMap()->getLocalNumElements());
1813 for (
size_t i = 0; i < (size_t)dirichletRows.
size(); i++) {
1814 if (dirichletRows[i]) {
1820 for (
size_t j = 0; j < (size_t)indices.
size(); j++)
1821 valuesNC[j] = replaceWith;
1826 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1831 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.
size()) == X->getMap()->getLocalNumElements());
1832 for (
size_t i = 0; i < (size_t)dirichletRows.
size(); i++) {
1833 if (dirichletRows[i]) {
1834 for (
size_t j = 0; j < X->getNumVectors(); j++)
1835 X->replaceLocalValue(i, j, replaceWith);
1840 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1843 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
1845 using ATS = Kokkos::ArithTraits<Scalar>;
1846 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1848 typename ATS::val_type impl_replaceWith = replaceWith;
1850 auto localMatrix = A->getLocalMatrixDevice();
1853 Kokkos::parallel_for(
1854 "MueLu:Utils::ZeroDirichletRows", range_type(0, numRows),
1856 if (dirichletRows(row)) {
1857 auto rowView = localMatrix.row(row);
1858 auto length = rowView.length;
1859 for (decltype(length) colID = 0; colID < length; colID++)
1860 rowView.value(colID) = impl_replaceWith;
1865 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1868 const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows,
1870 using ATS = Kokkos::ArithTraits<Scalar>;
1871 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1873 typename ATS::val_type impl_replaceWith = replaceWith;
1875 auto myCols = X->getLocalViewDevice(Xpetra::Access::ReadWrite);
1876 size_t numVecs = X->getNumVectors();
1877 Kokkos::parallel_for(
1878 "MueLu:Utils::ZeroDirichletRows_MV", range_type(0, dirichletRows.size()),
1879 KOKKOS_LAMBDA(
const size_t i) {
1880 if (dirichletRows(i)) {
1881 for (
size_t j = 0; j < numVecs; j++)
1882 myCols(i, j) = impl_replaceWith;
1887 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1892 TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.
size()) == A->getColMap()->getLocalNumElements());
1893 for (
size_t i = 0; i < A->getLocalNumRows(); i++) {
1896 A->getLocalRowView(i, indices, values);
1899 for (
size_t j = 0; j < static_cast<size_t>(indices.
size()); j++)
1900 if (dirichletCols[indices[j]])
1901 valuesNC[j] = replaceWith;
1905 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1908 const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols,
1910 using ATS = Kokkos::ArithTraits<Scalar>;
1911 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1913 typename ATS::val_type impl_replaceWith = replaceWith;
1915 auto localMatrix = A->getLocalMatrixDevice();
1918 Kokkos::parallel_for(
1919 "MueLu:Utils::ZeroDirichletCols", range_type(0, numRows),
1921 auto rowView = localMatrix.row(row);
1922 auto length = rowView.length;
1923 for (decltype(length) colID = 0; colID < length; colID++)
1924 if (dirichletCols(rowView.colidx(colID))) {
1925 rowView.value(colID) = impl_replaceWith;
1930 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1937 throw std::runtime_error(
"UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
1940 bool has_import = !importer.
is_null();
1943 std::vector<LocalOrdinal> dirichletRows;
1944 FindDirichletRows(A, dirichletRows);
1947 printf(
"[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
1948 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++)
1949 printf(
"%d ",dirichletRows[i]);
1962 for (
size_t i = 0; i < (size_t)dirichletRows.size(); i++) {
1963 dr[dirichletRows[i]] = 1;
1964 if (!has_import) dc[dirichletRows[i]] = 1;
1969 isDirichletCol->doImport(*isDirichletRow, *importer, Xpetra::CombineMode::ADD);
1972 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1976 using ISC =
typename Kokkos::ArithTraits<Scalar>::val_type;
1977 using range_type = Kokkos::RangePolicy<LocalOrdinal, typename Node::execution_space>;
1978 using local_matrix_type =
typename CrsMatrix::local_matrix_type;
1979 using values_type =
typename local_matrix_type::values_type;
1981 const ISC ONE = Kokkos::ArithTraits<ISC>::one();
1982 const ISC
ZERO = Kokkos::ArithTraits<ISC>::zero();
1985 auto localMatrix = original->getLocalMatrixDevice();
1987 values_type new_values(
"values", localMatrix.nnz());
1989 Kokkos::parallel_for(
1990 "ReplaceNonZerosWithOnes", range_type(0, localMatrix.nnz()), KOKKOS_LAMBDA(
const size_t i) {
1991 if (localMatrix.values(i) != ZERO)
1992 new_values(i) = ONE;
1994 new_values(i) = ZERO;
2000 NewMatrix->fillComplete(original->getDomainMap(), original->getRangeMap());
2004 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2015 if (!stridedMap.
is_null()) fullMap = stridedMap->getMap();
2018 const size_t numSubMaps = sourceBlockedMap.
getNumMaps();
2020 throw std::runtime_error(
"GenerateBlockedTargetMap(): Map compatibility error");
2025 for (
size_t i = 0; i < numSubMaps; i++) {
2028 for (
size_t j = 0; j < map->getLocalNumElements(); j++) {
2029 LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
2030 block_ids->replaceLocalValue(jj, (
int)i);
2037 new_block_ids->doImport(*block_ids, Importer, Xpetra::CombineMode::ADD);
2043 for (
size_t i = 0; i < targetMap->getLocalNumElements(); i++) {
2044 elementsInSubMap[data[i]].
push_back(targetMap->getGlobalElement(i));
2048 std::vector<RCP<const Map>> subMaps(numSubMaps);
2049 for (
size_t i = 0; i < numSubMaps; i++) {
2054 return rcp(
new BlockedMap(targetMap, subMaps));
2057 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2063 return tpColMap.isLocallyFitted(tpRowMap);
2069 const size_t numElements = rowElements.
size();
2071 if (
size_t(colElements.
size()) < numElements)
2074 bool goodMap =
true;
2075 for (
size_t i = 0; i < numElements; i++)
2076 if (rowElements[i] != colElements[i]) {
2084 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2089 using local_graph_type =
typename local_matrix_type::staticcrsgraph_type;
2090 using lno_nnz_view_t =
typename local_graph_type::entries_type::non_const_type;
2091 using device =
typename local_graph_type::device_type;
2092 using execution_space =
typename local_matrix_type::execution_space;
2093 using ordinal_type =
typename local_matrix_type::ordinal_type;
2095 local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
2097 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);
2103 auto view1D = Kokkos::subview(retval->getLocalViewDevice(Xpetra::Access::ReadWrite), Kokkos::ALL(), 0);
2104 Kokkos::parallel_for(
2105 "Utilities::ReverseCuthillMcKee",
2106 Kokkos::RangePolicy<ordinal_type, execution_space>(0, localGraph.numRows()),
2107 KOKKOS_LAMBDA(
const ordinal_type rowIdx) {
2108 view1D(rcmOrder(rowIdx)) = rowIdx;
2113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2118 using local_graph_type =
typename local_matrix_type::staticcrsgraph_type;
2119 using lno_nnz_view_t =
typename local_graph_type::entries_type::non_const_type;
2120 using device =
typename local_graph_type::device_type;
2121 using execution_space =
typename local_matrix_type::execution_space;
2122 using ordinal_type =
typename local_matrix_type::ordinal_type;
2124 local_graph_type localGraph = Op.getLocalMatrixDevice().graph;
2127 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);
2133 auto view1D = Kokkos::subview(retval->getLocalViewDevice(Xpetra::Access::ReadWrite), Kokkos::ALL(), 0);
2135 Kokkos::parallel_for(
2136 "Utilities::ReverseCuthillMcKee",
2137 Kokkos::RangePolicy<ordinal_type, execution_space>(0, numRows),
2138 KOKKOS_LAMBDA(
const ordinal_type rowIdx) {
2139 view1D(rcmOrder(numRows - 1 - rowIdx)) = rowIdx;
2146 #define MUELU_UTILITIESBASE_SHORT
2147 #endif // MUELU_UTILITIESBASE_DEF_HPP
static void FindNonZeros(const Teuchos::ArrayRCP< const Scalar > vals, Teuchos::ArrayRCP< bool > nonzeros)
Find non-zero values in an ArrayRCP Compares the value to 2 * machine epsilon.
static RCP< Export< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
static RCP< Vector > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const std::vector< LocalOrdinal > &dirichletRows)
static Teuchos::RCP< Vector > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i=k}(-a_ik), for each for i in the matrix.
static void ZeroDirichletCols(Teuchos::RCP< Matrix > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomain)
Detects Dirichlet columns & domains from a list of Dirichlet rows.
#define MueLu_sumAll(rcpComm, in, out)
Array< T > & append(const T &x)
virtual int getSize() const =0
MueLu::DefaultLocalOrdinal LocalOrdinal
void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, memory_space > &dirichletRows)
virtual LO getBlockSize() const override
static bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Matrix &Op)
static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node >> &isDirichletRow, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node >> &isDirichletCol)
virtual int getRank() const =0
static void ApplyRowSumCriterionHost(const Matrix &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, Kokkos::HostSpace > &dirichletRows)
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
static magnitudeType eps()
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)
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Apply Rowsum Criterion.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static RCP< Vector > GetMatrixOverlappedDeletedRowsum(const Matrix &A)
Extract Overlapped Matrix Deleted Rowsum.
virtual const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
static RCP< Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > GetThresholdedGraph(const RCP< Matrix > &A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a graph.
virtual void resumeFill(const RCP< ParameterList > ¶ms=null)=0
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > ReplaceNonZerosWithOnes(const RCP< Matrix > &original)
Creates a copy of a matrix where the non-zero entries are replaced by ones.
size_t getLocalNumRows() const override
Exception throws to report incompatible objects (like maps).
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Matrix &Op)
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getTargetMap() const =0
static RCP< Time > getNewTimer(const std::string &name)
void resize(const size_type n, const T &val=T())
virtual void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > ¶ms=null)=0
virtual void getLocalDiagCopy(Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &diag) const =0
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static Teuchos::ArrayRCP< const bool > DetectDirichletCols(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
virtual void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const =0
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
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)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static RCP< CrsMatrixWrap > GetThresholdedMatrix(const RCP< Matrix > &Ain, const Magnitude threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a matrix.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
virtual bool isFillComplete() const =0
static RCP< Xpetra::Vector< Magnitude, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedAbsDeletedRowsum(const Matrix &A)
size_t getNumMaps() const
virtual void replaceLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)=0
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal_arcp(const Matrix &A)
Extract Matrix Diagonal.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
Teuchos::ScalarTraits< Scalar >::magnitudeType Magnitude
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
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.
static void seedrandom(unsigned int s)
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows (extended version)
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
const RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getMap(size_t i, bool bThyraMode=false) const
virtual UnderlyingLib lib() const
static magnitudeType magnitude(T a)
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
void push_back(const value_type &x)
virtual RCP< const CrsGraph > getCrsGraph() const =0
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static void EnforceInitialCondition(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &InitialGuess, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), const bool count_twos_as_dirichlet=false)
Detect Dirichlet rows and copy values from RHS multivector to InitialGuess for Dirichlet rows...
impl_scalar_type_dualview::t_dev::const_type getValuesDevice(const LO &lclRow) const
virtual Teuchos::RCP< const Map > getRangeMap() const =0
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
Exception throws to report errors in the internal logical of the program.
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)
Frobenius inner product of two matrices.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getSourceMap() const =0
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
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.
static Kokkos::View< bool *, typename Kokkos::HostSpace > DetectDirichletRows_kokkos_host(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
virtual size_t getNumVectors() const =0
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
virtual UnderlyingLib lib() const =0
static Kokkos::View< bool *, typename NO::device_type::memory_space > DetectDirichletRows_kokkos(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
static RCP< Matrix > Crs2Op(RCP< CrsMatrix > Op)
virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const =0
virtual Teuchos::RCP< const Map > getDomainMap() const =0
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)