10 #ifndef MUELU_SAPFACTORY_DEF_HPP
11 #define MUELU_SAPFACTORY_DEF_HPP
13 #if KOKKOS_VERSION >= 40799
14 #include "KokkosKernels_ArithTraits.hpp"
16 #include "Kokkos_ArithTraits.hpp"
27 #include "MueLu_PerfUtils.hpp"
28 #include "MueLu_TentativePFactory.hpp"
29 #include "MueLu_Utilities.hpp"
35 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
38 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
41 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
45 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
58 #undef SET_VALID_ENTRY
60 validParamList->
set<
RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
65 norecurse.disableRecursiveValidation();
66 validParamList->
set<
ParameterList>(
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
68 return validParamList;
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
73 Input(fineLevel,
"A");
78 if (initialPFact == Teuchos::null) {
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
86 return BuildP(fineLevel, coarseLevel);
89 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 const std::string prefix =
"MueLu::SaPFactory(" + levelIDs +
"): ";
104 if (initialPFact == Teuchos::null) {
115 if (Ptent == Teuchos::null) {
116 finalP = Teuchos::null;
117 Set(coarseLevel,
"P", finalP);
121 if (restrictionMode_) {
130 if (pL.
isSublist(
"matrixmatrix: kernel params"))
134 if (coarseLevel.
IsAvailable(
"AP reuse data",
this)) {
135 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
143 APparams->
set(
"compute global constants: temporaries", APparams->
get(
"compute global constants: temporaries",
false));
144 APparams->
set(
"compute global constants", APparams->
get(
"compute global constants",
false));
146 const SC dampingFactor = as<SC>(pL.
get<
double>(
"sa: damping factor"));
147 const LO maxEigenIterations = as<LO>(pL.
get<
int>(
"sa: eigenvalue estimate num iterations"));
148 const bool estimateMaxEigen = pL.
get<
bool>(
"sa: calculate eigenvalue estimate");
149 const bool useAbsValueRowSum = pL.
get<
bool>(
"sa: use rowsumabs diagonal scaling");
150 const bool doQRStep = pL.
get<
bool>(
"tentative: calculate qr");
151 const bool enforceConstraints = pL.
get<
bool>(
"sa: enforce constraints");
152 const MT userDefinedMaxEigen = as<MT>(pL.
get<
double>(
"sa: max eigenvalue"));
153 double dTol = pL.
get<
double>(
"sa: rowsumabs diagonal replacement tolerance");
154 const MT diagonalReplacementTolerance = (dTol == as<double>(-1) ?
Teuchos::ScalarTraits<MT>::eps() * 100 : as<MT>(pL.
get<
double>(
"sa: rowsumabs diagonal replacement tolerance")));
155 const SC diagonalReplacementValue = as<SC>(pL.
get<
double>(
"sa: rowsumabs diagonal replacement value"));
156 const bool replaceSingleEntryRowWithZero = pL.
get<
bool>(
"sa: rowsumabs replace single entry row with zero");
157 const bool useAutomaticDiagTol = pL.
get<
bool>(
"sa: rowsumabs use automatic diagonal tolerance");
161 "MueLu::TentativePFactory::MakeTentative: cannot use 'enforce constraints' and 'calculate qr' at the same time");
166 if (userDefinedMaxEigen == -1.) {
168 lambdaMax = A->GetMaxEigenvalueEstimate();
170 GetOStream(
Statistics1) <<
"Calculating max eigenvalue estimate now (max iters = " << maxEigenIterations << ((useAbsValueRowSum) ?
", use rowSumAbs diagonal)" :
", use point diagonal)") << std::endl;
171 Coordinate stopTol = 1e-4;
172 if (useAbsValueRowSum) {
173 const bool returnReciprocal =
true;
175 diagonalReplacementTolerance,
176 diagonalReplacementValue,
177 replaceSingleEntryRowWithZero,
178 useAutomaticDiagTol);
180 "SaPFactory: eigenvalue estimate: diagonal reciprocal is null.");
184 A->SetMaxEigenvalueEstimate(lambdaMax);
186 GetOStream(
Statistics1) <<
"Using cached max eigenvalue estimate" << std::endl;
189 lambdaMax = userDefinedMaxEigen;
190 A->SetMaxEigenvalueEstimate(lambdaMax);
192 GetOStream(
Statistics1) <<
"Prolongator damping factor = " << dampingFactor / lambdaMax <<
" (" << dampingFactor <<
" / " << lambdaMax <<
")" << std::endl;
196 if (!useAbsValueRowSum)
198 else if (invDiag == Teuchos::null) {
199 GetOStream(
Runtime0) <<
"Using rowsumabs diagonal" << std::endl;
200 const bool returnReciprocal =
true;
202 diagonalReplacementTolerance,
203 diagonalReplacementValue,
204 replaceSingleEntryRowWithZero,
205 useAutomaticDiagTol);
209 SC omega = dampingFactor / lambdaMax;
215 finalP =
Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(
Statistics2), std::string(
"MueLu::SaP-") +
toString(coarseLevel.
GetLevelID()), APparams);
216 if (enforceConstraints) {
217 if (!pL.
get<
bool>(
"use kokkos refactor")) {
218 if (A->GetFixedBlockSize() == 1)
219 optimalSatisfyPConstraintsForScalarPDEsNonKokkos(finalP);
221 SatisfyPConstraintsNonKokkos(A, finalP);
226 SatisfyPConstraints(A, finalP);
238 if (!restrictionMode_) {
241 std::ostringstream oss;
243 finalP->setObjectLabel(oss.str());
245 Set(coarseLevel,
"P", finalP);
247 APparams->
set(
"graph", finalP);
248 Set(coarseLevel,
"AP reuse data", APparams);
251 if (Ptent->IsView(
"stridedMaps"))
252 finalP->CreateView(
"stridedMaps", Ptent);
260 std::ostringstream oss;
262 R->setObjectLabel(oss.str());
266 Set(coarseLevel,
"R", R);
269 if (Ptent->IsView(
"stridedMaps"))
270 R->CreateView(
"stridedMaps", Ptent,
true );
275 params->
set(
"printLoadBalancingInfo",
true);
276 params->
set(
"printCommInfo",
true);
298 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
302 LO nPDEs = A->GetFixedBlockSize();
306 for (
size_t k = 0; k < (size_t)nPDEs; k++) ConstraintViolationSum[k] = zero;
307 for (
size_t k = 0; k < (size_t)nPDEs; k++) Rsum[k] = zero;
308 for (
size_t k = 0; k < (size_t)nPDEs; k++) nPositive[k] = 0;
310 for (
size_t i = 0; i < as<size_t>(P->getRowMap()->getLocalNumElements()); i++) {
314 P->getLocalRowView((
LO)i, indices, vals1);
315 size_t nnz = indices.
size();
316 if (nnz == 0)
continue;
320 bool checkRow =
true;
325 for (
LO j = 0; j < indices.
size(); j++) {
326 Rsum[j % nPDEs] += vals[j];
328 ConstraintViolationSum[j % nPDEs] += vals[j];
332 (nPositive[j % nPDEs])++;
335 ConstraintViolationSum[j % nPDEs] += (vals[j] - one);
345 for (
size_t k = 0; k < (size_t)nPDEs; k++) {
347 ConstraintViolationSum[k] += (-Rsum[k]);
349 ConstraintViolationSum[k] += (one - Rsum[k]);
354 for (
size_t k = 0; k < (size_t)nPDEs; k++) {
360 for (
LO j = 0; j < indices.
size(); j++) {
362 vals[j] += (ConstraintViolationSum[j % nPDEs] / (as<Scalar>(nPositive[j % nPDEs])));
365 for (
size_t k = 0; k < (size_t)nPDEs; k++) ConstraintViolationSum[k] = zero;
367 for (
size_t k = 0; k < (size_t)nPDEs; k++) Rsum[k] = zero;
368 for (
size_t k = 0; k < (size_t)nPDEs; k++) nPositive[k] = 0;
373 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
382 for (
size_t i = 0; i < as<size_t>(P->getRowMap()->getLocalNumElements()); i++) {
387 size_t nnz = indices.
size();
391 for (
size_t j = 0; j < nnz; j++) rsumTarget += vals[j];
393 if (nnz > as<size_t>(maxEntriesPerRow)) {
394 maxEntriesPerRow = nnz * 3;
395 scalarData.
resize(3 * maxEntriesPerRow);
397 hasFeasible = constrainRow(vals.
getRawPtr(), as<LocalOrdinal>(nnz), zero, one, rsumTarget, vals.
getRawPtr(), scalarData.
getRawPtr());
400 for (
size_t j = 0; j < nnz; j++) vals[j] = one / as<Scalar>(nnz);
407 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
526 Scalar notFlippedLeftBound, notFlippedRghtBound, aBigNumber, *origSorted;
527 Scalar rowSumDeviation, temp, *fixedSorted, delta;
528 Scalar closestToLeftBoundDist, closestToRghtBoundDist;
533 notFlippedLeftBound = leftBound;
534 notFlippedRghtBound = rghtBound;
538 hasFeasibleSol =
true;
540 hasFeasibleSol =
false;
541 return hasFeasibleSol;
553 origSorted = &scalarData[0];
554 fixedSorted = &(scalarData[nEntries]);
557 for (
LocalOrdinal i = 0; i < nEntries; i++) inds[i] = i;
558 for (
LocalOrdinal i = 0; i < nEntries; i++) origSorted[i] = orig[i];
561 std::sort(inds, inds + nEntries,
564 for (
LocalOrdinal i = 0; i < nEntries; i++) origSorted[i] = orig[inds[i]];
566 closestToLeftBound = 0;
570 closestToRghtBound = closestToLeftBound;
576 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
577 if (closestToRghtBound == nEntries)
578 closestToRghtBoundDist = aBigNumber;
580 closestToRghtBoundDist = origSorted[closestToRghtBound] - rghtBound;
585 rowSumDeviation = leftBound * as<Scalar>(closestToLeftBound) + as<Scalar>((nEntries - closestToRghtBound)) * rghtBound - rsumTarget;
586 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) rowSumDeviation += origSorted[i];
594 leftBound = -rghtBound;
599 if ((nEntries % 2) == 1) origSorted[(nEntries / 2)] = -origSorted[(nEntries / 2)];
601 temp = origSorted[i];
602 origSorted[i] = -origSorted[nEntries - 1 - i];
603 origSorted[nEntries - i - 1] = -temp;
609 closestToLeftBound = nEntries - closestToRghtBound;
610 closestToRghtBound = nEntries - itemp;
611 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
612 if (closestToRghtBound == nEntries)
613 closestToRghtBoundDist = aBigNumber;
615 closestToRghtBoundDist = origSorted[closestToRghtBound] - rghtBound;
617 rowSumDeviation = -rowSumDeviation;
622 for (
LocalOrdinal i = 0; i < closestToLeftBound; i++) fixedSorted[i] = leftBound;
623 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) fixedSorted[i] = origSorted[i];
624 for (
LocalOrdinal i = closestToRghtBound; i < nEntries; i++) fixedSorted[i] = rghtBound;
627 if (closestToRghtBound != closestToLeftBound)
628 delta = rowSumDeviation / as<Scalar>(closestToRghtBound - closestToLeftBound);
635 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) fixedSorted[i] = origSorted[i] - delta;
637 rowSumDeviation = rowSumDeviation - closestToLeftBoundDist;
638 fixedSorted[closestToLeftBound] = leftBound;
639 closestToLeftBound++;
640 if (closestToLeftBound < nEntries)
641 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
643 closestToLeftBoundDist = aBigNumber;
648 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) fixedSorted[i] = origSorted[i] - delta;
650 rowSumDeviation = rowSumDeviation + closestToRghtBoundDist;
652 fixedSorted[closestToRghtBound] = origSorted[closestToRghtBound];
653 closestToRghtBound++;
655 if (closestToRghtBound >= nEntries)
656 closestToRghtBoundDist = aBigNumber;
658 closestToRghtBoundDist = origSorted[closestToRghtBound] - rghtBound;
666 if ((nEntries % 2) == 1) fixedSorted[(nEntries / 2)] = -fixedSorted[(nEntries / 2)];
668 temp = fixedSorted[i];
669 fixedSorted[i] = -fixedSorted[nEntries - 1 - i];
670 fixedSorted[nEntries - i - 1] = -temp;
673 for (
LocalOrdinal i = 0; i < nEntries; i++) fixedUnsorted[inds[i]] = fixedSorted[i];
677 bool lowerViolation =
false;
678 bool upperViolation =
false;
679 bool sumViolation =
false;
683 if (TST::real(fixedUnsorted[i]) < TST::real(notFlippedLeftBound)) lowerViolation =
true;
684 if (TST::real(fixedUnsorted[i]) > TST::real(notFlippedRghtBound)) upperViolation =
true;
685 temp += fixedUnsorted[i];
687 SC tol = as<Scalar>(std::max(1.0e-8, as<double>(100 * TST::eps())));
688 if (TST::magnitude(temp - rsumTarget) > TST::magnitude(tol * rsumTarget)) sumViolation =
true;
694 return hasFeasibleSol;
697 template <
typename local_matrix_type>
699 using Scalar =
typename local_matrix_type::non_const_value_type;
701 using LO =
typename local_matrix_type::non_const_ordinal_type;
702 using Device =
typename local_matrix_type::device_type;
703 #if KOKKOS_VERSION >= 40799
704 using KAT = KokkosKernels::ArithTraits<SC>;
706 using KAT = Kokkos::ArithTraits<SC>;
713 Kokkos::View<SC**, Device>
Rsum;
720 Rsum = Kokkos::View<SC**, Device>(
"Rsum", localP_.numRows(),
nPDEs);
721 nPositive = Kokkos::View<size_t**, Device>(
"nPositive", localP_.numRows(),
nPDEs);
724 KOKKOS_INLINE_FUNCTION
726 auto rowPtr =
localP.graph.row_map;
727 auto values =
localP.values;
729 bool checkRow =
true;
731 if (rowPtr(rowIdx + 1) == rowPtr(rowIdx)) checkRow =
false;
736 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) {
737 Rsum(rowIdx, entryIdx %
nPDEs) += values(entryIdx);
738 if (KAT::real(values(entryIdx)) < KAT::real(
zero)) {
740 values(entryIdx) =
zero;
742 if (KAT::real(values(entryIdx)) != KAT::real(
zero))
745 if (KAT::real(values(entryIdx)) > KAT::real(1.00001)) {
747 values(entryIdx) =
one;
756 for (
size_t k = 0; k < (size_t)
nPDEs; k++) {
757 if (KAT::real(
Rsum(rowIdx, k)) < KAT::magnitude(
zero)) {
759 }
else if (KAT::real(
Rsum(rowIdx, k)) > KAT::magnitude(1.00001)) {
765 for (
size_t k = 0; k < (size_t)
nPDEs; k++) {
771 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) {
772 if (KAT::real(values(entryIdx)) > KAT::real(
zero)) {
773 values(entryIdx) = values(entryIdx) +
779 for (
size_t k = 0; k < (size_t)
nPDEs; k++)
Rsum(rowIdx, k) =
zero;
780 for (
size_t k = 0; k < (size_t)
nPDEs; k++)
nPositive(rowIdx, k) = 0;
785 template <
typename local_matrix_type>
787 using Scalar =
typename local_matrix_type::non_const_value_type;
789 using LO =
typename local_matrix_type::non_const_ordinal_type;
790 using Device =
typename local_matrix_type::device_type;
791 #if KOKKOS_VERSION >= 40799
792 using KAT = KokkosKernels::ArithTraits<SC>;
794 using KAT = Kokkos::ArithTraits<SC>;
802 Kokkos::View<LO**, Device>
inds;
807 origSorted = Kokkos::View<SC**, Device>(
"origSorted", localP_.numRows(), maxRowEntries_);
808 fixedSorted = Kokkos::View<SC**, Device>(
"fixedSorted", localP_.numRows(), maxRowEntries_);
809 inds = Kokkos::View<LO**, Device>(
"inds", localP_.numRows(), maxRowEntries_);
812 KOKKOS_INLINE_FUNCTION
814 auto rowPtr =
localP.graph.row_map;
815 auto values =
localP.values;
817 auto nnz = rowPtr(rowIdx + 1) - rowPtr(rowIdx);
821 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) rsumTarget += values(entryIdx);
825 SC rowSumDeviation, temp, delta;
826 SC closestToLeftBoundDist, closestToRghtBoundDist;
827 LO closestToLeftBound, closestToRghtBound;
832 if ((KAT::real(rsumTarget) >= KAT::real(leftBound * (static_cast<SC>(nnz)))) &&
833 (KAT::real(rsumTarget) <= KAT::real(rghtBound * (static_cast<SC>(nnz))))) {
838 aBigNumber = KAT::zero();
839 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) {
840 if (KAT::magnitude(values(entryIdx)) > KAT::magnitude(aBigNumber))
841 aBigNumber = KAT::magnitude(values(entryIdx));
843 aBigNumber = aBigNumber + (KAT::magnitude(leftBound) + KAT::magnitude(rghtBound)) * (100 *
one);
846 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) {
848 inds(rowIdx, ind) = ind;
852 auto sortVals = Kokkos::subview(
origSorted, rowIdx, Kokkos::ALL());
853 auto sortInds = Kokkos::subview(
inds, rowIdx, Kokkos::make_pair(0, ind));
857 for (
LO i = 1; i < static_cast<LO>(nnz); ++i) {
861 if (KAT::real(sortVals(sortInds(i))) < KAT::real(sortVals(sortInds(0)))) {
862 for (; j > 0; --j) sortInds(j) = sortInds(j - 1);
866 for (; KAT::real(sortVals(ind)) < KAT::real(sortVals(sortInds(j - 1))); --j) sortInds(j) = sortInds(j - 1);
872 for (
LO i = 0; i < static_cast<LO>(nnz); i++)
origSorted(rowIdx, i) = values(rowPtr(rowIdx) +
inds(rowIdx, i));
874 closestToLeftBound = 0;
875 while ((closestToLeftBound < static_cast<LO>(nnz)) && (KAT::real(
origSorted(rowIdx, closestToLeftBound)) <= KAT::real(leftBound))) closestToLeftBound++;
878 closestToRghtBound = closestToLeftBound;
879 while ((closestToRghtBound < static_cast<LO>(nnz)) && (KAT::real(
origSorted(rowIdx, closestToRghtBound)) <= KAT::real(rghtBound))) closestToRghtBound++;
884 closestToLeftBoundDist =
origSorted(rowIdx, closestToLeftBound) - leftBound;
885 if (closestToRghtBound == static_cast<LO>(nnz))
886 closestToRghtBoundDist = aBigNumber;
888 closestToRghtBoundDist =
origSorted(rowIdx, closestToRghtBound) - rghtBound;
893 rowSumDeviation = leftBound * (
static_cast<SC>(closestToLeftBound)) + (
static_cast<SC>(nnz - closestToRghtBound)) * rghtBound - rsumTarget;
894 for (
LO i = closestToLeftBound; i < closestToRghtBound; i++) rowSumDeviation +=
origSorted(rowIdx, i);
899 if (KAT::real(rowSumDeviation) < KAT::real(KAT::zero())) {
902 leftBound = -rghtBound;
909 for (
LO i = 0; i < static_cast<LO>(nnz / 2); i++) {
917 LO itemp = closestToLeftBound;
918 closestToLeftBound = nnz - closestToRghtBound;
919 closestToRghtBound = nnz - itemp;
920 closestToLeftBoundDist =
origSorted(rowIdx, closestToLeftBound) - leftBound;
921 if (closestToRghtBound == static_cast<LO>(nnz))
922 closestToRghtBoundDist = aBigNumber;
924 closestToRghtBoundDist =
origSorted(rowIdx, closestToRghtBound) - rghtBound;
926 rowSumDeviation = -rowSumDeviation;
931 for (
LO i = 0; i < closestToLeftBound; i++)
fixedSorted(rowIdx, i) = leftBound;
933 for (
LO i = closestToRghtBound; i < static_cast<LO>(nnz); i++)
fixedSorted(rowIdx, i) = rghtBound;
935 while ((KAT::magnitude(rowSumDeviation) > KAT::magnitude((
one * 1.e-10) * rsumTarget))) {
936 if (closestToRghtBound != closestToLeftBound)
937 delta = rowSumDeviation /
static_cast<SC>(closestToRghtBound - closestToLeftBound);
941 if (KAT::magnitude(closestToLeftBoundDist) <= KAT::magnitude(closestToRghtBoundDist)) {
942 if (KAT::magnitude(delta) <= KAT::magnitude(closestToLeftBoundDist)) {
943 rowSumDeviation =
zero;
944 for (
LO i = closestToLeftBound; i < closestToRghtBound; i++)
fixedSorted(rowIdx, i) =
origSorted(rowIdx, i) - delta;
946 rowSumDeviation = rowSumDeviation - closestToLeftBoundDist;
947 fixedSorted(rowIdx, closestToLeftBound) = leftBound;
948 closestToLeftBound++;
949 if (closestToLeftBound < static_cast<LO>(nnz))
950 closestToLeftBoundDist =
origSorted(rowIdx, closestToLeftBound) - leftBound;
952 closestToLeftBoundDist = aBigNumber;
955 if (KAT::magnitude(delta) <= KAT::magnitude(closestToRghtBoundDist)) {
957 for (
LO i = closestToLeftBound; i < closestToRghtBound; i++)
fixedSorted(rowIdx, i) =
origSorted(rowIdx, i) - delta;
959 rowSumDeviation = rowSumDeviation + closestToRghtBoundDist;
962 closestToRghtBound++;
964 if (closestToRghtBound >= static_cast<LO>(nnz))
965 closestToRghtBoundDist = aBigNumber;
967 closestToRghtBoundDist =
origSorted(rowIdx, closestToRghtBound) - rghtBound;
972 auto rowStart = rowPtr(rowIdx);
977 for (
LO i = 0; i < static_cast<LO>(nnz / 2); i++) {
984 for (
LO i = 0; i < static_cast<LO>(nnz); i++) values(rowStart +
inds(rowIdx, i)) =
fixedSorted(rowIdx, i);
988 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) values(entryIdx) =
one / (
static_cast<SC>(nnz));
995 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
997 using Device =
typename Matrix::local_matrix_type::device_type;
998 LO nPDEs = A->GetFixedBlockSize();
1000 using local_mat_type =
typename Matrix::local_matrix_type;
1002 Kokkos::parallel_for(
"enforce constraint", Kokkos::RangePolicy<typename Device::execution_space>(0, P->getRowMap()->getLocalNumElements()),
1007 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1009 using Device =
typename Matrix::local_matrix_type::device_type;
1012 using local_mat_type =
typename Matrix::local_matrix_type;
1014 Kokkos::parallel_for(
"enforce constraint", Kokkos::RangePolicy<typename Device::execution_space>(0, P->getLocalNumRows()),
1021 #endif // MUELU_SAPFACTORY_DEF_HPP
typename local_matrix_type::non_const_value_type Scalar
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
void optimalSatisfyPConstraintsForScalarPDEs(RCP< Matrix > &P) const
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
typename local_matrix_type::device_type Device
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.
bool constrainRow(Scalar *orig, LocalOrdinal nEntries, Scalar leftBound, Scalar rghtBound, Scalar rsumTarget, Scalar *fixedUnsorted, Scalar *scalarData) const
optimalSatisfyConstraintsForScalarPDEsKernel(LO nPDEs_, LO maxRowEntries_, local_matrix_type localP_)
T & get(const std::string &name, T def_value)
KOKKOS_INLINE_FUNCTION void operator()(const size_t rowIdx) const
SaPFactory()
Constructor. User can supply a factory for generating the tentative prolongator.
Kokkos::View< LO **, Device > inds
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static magnitudeType real(T a)
typename local_matrix_type::non_const_ordinal_type LO
constraintKernel(LO nPDEs_, local_matrix_type localP_)
One-liner description of what is happening.
void SatisfyPConstraintsNonKokkos(RCP< Matrix > A, RCP< Matrix > &P) const
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Kokkos::View< SC **, Device > fixedSorted
Kokkos::ArithTraits< SC > KAT
Print even more statistics.
Kokkos::View< size_t **, Device > nPositive
void resize(const size_type n, const T &val=T())
bool isParameter(const std::string &name) const
typename local_matrix_type::device_type Device
static RCP< Matrix > Jacobi(SC omega, const Vector &Dinv, const Matrix &A, const Matrix &B, RCP< Matrix > C_in, Teuchos::FancyOStream &fos, const std::string &label, RCP< ParameterList > ¶ms)
virtual ~SaPFactory()
Destructor.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
Kokkos::View< SC **, Device > origSorted
bool isSublist(const std::string &name) const
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
Kokkos::View< SC **, Device > ConstraintViolationSum
#define SET_VALID_ENTRY(name)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void SatisfyPConstraints(RCP< Matrix > A, RCP< Matrix > &P) const
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.
Kokkos::ArithTraits< SC > KAT
static magnitudeType magnitude(T a)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
typename local_matrix_type::non_const_value_type Scalar
typename local_matrix_type::non_const_ordinal_type LO
Kokkos::View< SC **, Device > Rsum
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
KOKKOS_INLINE_FUNCTION void operator()(const size_t rowIdx) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
void optimalSatisfyPConstraintsForScalarPDEsNonKokkos(RCP< Matrix > &P) const
int GetLevelID() const
Return level number.
Exception throws to report errors in the internal logical of the program.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
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.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.