46 #ifndef MUELU_SAPFACTORY_DEF_HPP
47 #define MUELU_SAPFACTORY_DEF_HPP
49 #include "Kokkos_ArithTraits.hpp"
59 #include "MueLu_PerfUtils.hpp"
60 #include "MueLu_TentativePFactory.hpp"
61 #include "MueLu_Utilities.hpp"
67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
84 #undef SET_VALID_ENTRY
86 validParamList->
set<
RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
91 norecurse.disableRecursiveValidation();
92 validParamList->
set<
ParameterList>(
"matrixmatrix: kernel params", norecurse,
"MatrixMatrix kernel parameters");
94 return validParamList;
97 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
99 Input(fineLevel,
"A");
104 if (initialPFact == Teuchos::null) {
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
112 return BuildP(fineLevel, coarseLevel);
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
121 const std::string prefix =
"MueLu::SaPFactory(" + levelIDs +
"): ";
130 if (initialPFact == Teuchos::null) {
141 if (Ptent == Teuchos::null) {
142 finalP = Teuchos::null;
143 Set(coarseLevel,
"P", finalP);
147 if (restrictionMode_) {
156 if (pL.
isSublist(
"matrixmatrix: kernel params"))
160 if (coarseLevel.
IsAvailable(
"AP reuse data",
this)) {
161 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
169 APparams->
set(
"compute global constants: temporaries", APparams->
get(
"compute global constants: temporaries",
false));
170 APparams->
set(
"compute global constants", APparams->
get(
"compute global constants",
false));
172 const SC dampingFactor = as<SC>(pL.
get<
double>(
"sa: damping factor"));
173 const LO maxEigenIterations = as<LO>(pL.
get<
int>(
"sa: eigenvalue estimate num iterations"));
174 const bool estimateMaxEigen = pL.
get<
bool>(
"sa: calculate eigenvalue estimate");
175 const bool useAbsValueRowSum = pL.
get<
bool>(
"sa: use rowsumabs diagonal scaling");
176 const bool doQRStep = pL.
get<
bool>(
"tentative: calculate qr");
177 const bool enforceConstraints = pL.
get<
bool>(
"sa: enforce constraints");
178 const MT userDefinedMaxEigen = as<MT>(pL.
get<
double>(
"sa: max eigenvalue"));
179 double dTol = pL.
get<
double>(
"sa: rowsumabs diagonal replacement tolerance");
180 const MT diagonalReplacementTolerance = (dTol == as<double>(-1) ?
Teuchos::ScalarTraits<MT>::eps() * 100 : as<MT>(pL.
get<
double>(
"sa: rowsumabs diagonal replacement tolerance")));
181 const SC diagonalReplacementValue = as<SC>(pL.
get<
double>(
"sa: rowsumabs diagonal replacement value"));
182 const bool replaceSingleEntryRowWithZero = pL.
get<
bool>(
"sa: rowsumabs replace single entry row with zero");
183 const bool useAutomaticDiagTol = pL.
get<
bool>(
"sa: rowsumabs use automatic diagonal tolerance");
187 "MueLu::TentativePFactory::MakeTentative: cannot use 'enforce constraints' and 'calculate qr' at the same time");
192 if (userDefinedMaxEigen == -1.) {
194 lambdaMax = A->GetMaxEigenvalueEstimate();
196 GetOStream(
Statistics1) <<
"Calculating max eigenvalue estimate now (max iters = " << maxEigenIterations << ((useAbsValueRowSum) ?
", use rowSumAbs diagonal)" :
", use point diagonal)") << std::endl;
197 Coordinate stopTol = 1e-4;
198 if (useAbsValueRowSum) {
199 const bool returnReciprocal =
true;
201 diagonalReplacementTolerance,
202 diagonalReplacementValue,
203 replaceSingleEntryRowWithZero,
204 useAutomaticDiagTol);
206 "SaPFactory: eigenvalue estimate: diagonal reciprocal is null.");
210 A->SetMaxEigenvalueEstimate(lambdaMax);
212 GetOStream(
Statistics1) <<
"Using cached max eigenvalue estimate" << std::endl;
215 lambdaMax = userDefinedMaxEigen;
216 A->SetMaxEigenvalueEstimate(lambdaMax);
218 GetOStream(
Statistics1) <<
"Prolongator damping factor = " << dampingFactor / lambdaMax <<
" (" << dampingFactor <<
" / " << lambdaMax <<
")" << std::endl;
222 if (!useAbsValueRowSum)
224 else if (invDiag == Teuchos::null) {
225 GetOStream(
Runtime0) <<
"Using rowsumabs diagonal" << std::endl;
226 const bool returnReciprocal =
true;
228 diagonalReplacementTolerance,
229 diagonalReplacementValue,
230 replaceSingleEntryRowWithZero,
231 useAutomaticDiagTol);
235 SC omega = dampingFactor / lambdaMax;
241 finalP =
Xpetra::IteratorOps<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP, GetOStream(
Statistics2), std::string(
"MueLu::SaP-") +
toString(coarseLevel.
GetLevelID()), APparams);
242 if (enforceConstraints) {
243 if (!pL.
get<
bool>(
"use kokkos refactor")) {
244 if (A->GetFixedBlockSize() == 1)
245 optimalSatisfyPConstraintsForScalarPDEsNonKokkos(finalP);
247 SatisfyPConstraintsNonKokkos(A, finalP);
252 SatisfyPConstraints(A, finalP);
264 if (!restrictionMode_) {
267 std::ostringstream oss;
269 finalP->setObjectLabel(oss.str());
271 Set(coarseLevel,
"P", finalP);
273 APparams->
set(
"graph", finalP);
274 Set(coarseLevel,
"AP reuse data", APparams);
277 if (Ptent->IsView(
"stridedMaps"))
278 finalP->CreateView(
"stridedMaps", Ptent);
286 std::ostringstream oss;
288 R->setObjectLabel(oss.str());
292 Set(coarseLevel,
"R", R);
295 if (Ptent->IsView(
"stridedMaps"))
296 R->CreateView(
"stridedMaps", Ptent,
true );
301 params->
set(
"printLoadBalancingInfo",
true);
302 params->
set(
"printCommInfo",
true);
324 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
328 LO nPDEs = A->GetFixedBlockSize();
332 for (
size_t k = 0; k < (size_t)nPDEs; k++) ConstraintViolationSum[k] = zero;
333 for (
size_t k = 0; k < (size_t)nPDEs; k++) Rsum[k] = zero;
334 for (
size_t k = 0; k < (size_t)nPDEs; k++) nPositive[k] = 0;
336 for (
size_t i = 0; i < as<size_t>(P->getRowMap()->getLocalNumElements()); i++) {
340 P->getLocalRowView((
LO)i, indices, vals1);
341 size_t nnz = indices.
size();
342 if (nnz == 0)
continue;
346 bool checkRow =
true;
351 for (
LO j = 0; j < indices.
size(); j++) {
352 Rsum[j % nPDEs] += vals[j];
354 ConstraintViolationSum[j % nPDEs] += vals[j];
358 (nPositive[j % nPDEs])++;
361 ConstraintViolationSum[j % nPDEs] += (vals[j] - one);
371 for (
size_t k = 0; k < (size_t)nPDEs; k++) {
373 ConstraintViolationSum[k] += (-Rsum[k]);
375 ConstraintViolationSum[k] += (one - Rsum[k]);
380 for (
size_t k = 0; k < (size_t)nPDEs; k++) {
386 for (
LO j = 0; j < indices.
size(); j++) {
388 vals[j] += (ConstraintViolationSum[j % nPDEs] / (as<Scalar>(nPositive[j % nPDEs])));
391 for (
size_t k = 0; k < (size_t)nPDEs; k++) ConstraintViolationSum[k] = zero;
393 for (
size_t k = 0; k < (size_t)nPDEs; k++) Rsum[k] = zero;
394 for (
size_t k = 0; k < (size_t)nPDEs; k++) nPositive[k] = 0;
399 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
408 for (
size_t i = 0; i < as<size_t>(P->getRowMap()->getLocalNumElements()); i++) {
413 size_t nnz = indices.
size();
417 for (
size_t j = 0; j < nnz; j++) rsumTarget += vals[j];
419 if (nnz > as<size_t>(maxEntriesPerRow)) {
420 maxEntriesPerRow = nnz * 3;
421 scalarData.
resize(3 * maxEntriesPerRow);
423 hasFeasible = constrainRow(vals.
getRawPtr(), as<LocalOrdinal>(nnz), zero, one, rsumTarget, vals.
getRawPtr(), scalarData.
getRawPtr());
426 for (
size_t j = 0; j < nnz; j++) vals[j] = one / as<Scalar>(nnz);
433 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
552 Scalar notFlippedLeftBound, notFlippedRghtBound, aBigNumber, *origSorted;
553 Scalar rowSumDeviation, temp, *fixedSorted, delta;
554 Scalar closestToLeftBoundDist, closestToRghtBoundDist;
559 notFlippedLeftBound = leftBound;
560 notFlippedRghtBound = rghtBound;
564 hasFeasibleSol =
true;
566 hasFeasibleSol =
false;
567 return hasFeasibleSol;
579 origSorted = &scalarData[0];
580 fixedSorted = &(scalarData[nEntries]);
583 for (
LocalOrdinal i = 0; i < nEntries; i++) inds[i] = i;
584 for (
LocalOrdinal i = 0; i < nEntries; i++) origSorted[i] = orig[i];
587 std::sort(inds, inds + nEntries,
590 for (
LocalOrdinal i = 0; i < nEntries; i++) origSorted[i] = orig[inds[i]];
592 closestToLeftBound = 0;
596 closestToRghtBound = closestToLeftBound;
602 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
603 if (closestToRghtBound == nEntries)
604 closestToRghtBoundDist = aBigNumber;
606 closestToRghtBoundDist = origSorted[closestToRghtBound] - rghtBound;
611 rowSumDeviation = leftBound * as<Scalar>(closestToLeftBound) + as<Scalar>((nEntries - closestToRghtBound)) * rghtBound - rsumTarget;
612 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) rowSumDeviation += origSorted[i];
620 leftBound = -rghtBound;
625 if ((nEntries % 2) == 1) origSorted[(nEntries / 2)] = -origSorted[(nEntries / 2)];
627 temp = origSorted[i];
628 origSorted[i] = -origSorted[nEntries - 1 - i];
629 origSorted[nEntries - i - 1] = -temp;
635 closestToLeftBound = nEntries - closestToRghtBound;
636 closestToRghtBound = nEntries - itemp;
637 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
638 if (closestToRghtBound == nEntries)
639 closestToRghtBoundDist = aBigNumber;
641 closestToRghtBoundDist = origSorted[closestToRghtBound] - rghtBound;
643 rowSumDeviation = -rowSumDeviation;
648 for (
LocalOrdinal i = 0; i < closestToLeftBound; i++) fixedSorted[i] = leftBound;
649 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) fixedSorted[i] = origSorted[i];
650 for (
LocalOrdinal i = closestToRghtBound; i < nEntries; i++) fixedSorted[i] = rghtBound;
653 if (closestToRghtBound != closestToLeftBound)
654 delta = rowSumDeviation / as<Scalar>(closestToRghtBound - closestToLeftBound);
661 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) fixedSorted[i] = origSorted[i] - delta;
663 rowSumDeviation = rowSumDeviation - closestToLeftBoundDist;
664 fixedSorted[closestToLeftBound] = leftBound;
665 closestToLeftBound++;
666 if (closestToLeftBound < nEntries)
667 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
669 closestToLeftBoundDist = aBigNumber;
674 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) fixedSorted[i] = origSorted[i] - delta;
676 rowSumDeviation = rowSumDeviation + closestToRghtBoundDist;
678 fixedSorted[closestToRghtBound] = origSorted[closestToRghtBound];
679 closestToRghtBound++;
681 if (closestToRghtBound >= nEntries)
682 closestToRghtBoundDist = aBigNumber;
684 closestToRghtBoundDist = origSorted[closestToRghtBound] - rghtBound;
692 if ((nEntries % 2) == 1) fixedSorted[(nEntries / 2)] = -fixedSorted[(nEntries / 2)];
694 temp = fixedSorted[i];
695 fixedSorted[i] = -fixedSorted[nEntries - 1 - i];
696 fixedSorted[nEntries - i - 1] = -temp;
699 for (
LocalOrdinal i = 0; i < nEntries; i++) fixedUnsorted[inds[i]] = fixedSorted[i];
703 bool lowerViolation =
false;
704 bool upperViolation =
false;
705 bool sumViolation =
false;
709 if (TST::real(fixedUnsorted[i]) < TST::real(notFlippedLeftBound)) lowerViolation =
true;
710 if (TST::real(fixedUnsorted[i]) > TST::real(notFlippedRghtBound)) upperViolation =
true;
711 temp += fixedUnsorted[i];
713 SC tol = as<Scalar>(std::max(1.0e-8, as<double>(100 * TST::eps())));
714 if (TST::magnitude(temp - rsumTarget) > TST::magnitude(tol * rsumTarget)) sumViolation =
true;
720 return hasFeasibleSol;
723 template <
typename local_matrix_type>
725 using Scalar =
typename local_matrix_type::non_const_value_type;
727 using LO =
typename local_matrix_type::non_const_ordinal_type;
728 using Device =
typename local_matrix_type::device_type;
729 using KAT = Kokkos::ArithTraits<SC>;
735 Kokkos::View<SC**, Device>
Rsum;
742 Rsum = Kokkos::View<SC**, Device>(
"Rsum", localP_.numRows(),
nPDEs);
743 nPositive = Kokkos::View<size_t**, Device>(
"nPositive", localP_.numRows(),
nPDEs);
746 KOKKOS_INLINE_FUNCTION
748 auto rowPtr =
localP.graph.row_map;
749 auto values =
localP.values;
751 bool checkRow =
true;
753 if (rowPtr(rowIdx + 1) == rowPtr(rowIdx)) checkRow =
false;
758 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) {
759 Rsum(rowIdx, entryIdx %
nPDEs) += values(entryIdx);
760 if (KAT::real(values(entryIdx)) < KAT::real(
zero)) {
762 values(entryIdx) =
zero;
764 if (KAT::real(values(entryIdx)) != KAT::real(
zero))
767 if (KAT::real(values(entryIdx)) > KAT::real(1.00001)) {
769 values(entryIdx) =
one;
778 for (
size_t k = 0; k < (size_t)
nPDEs; k++) {
779 if (KAT::real(
Rsum(rowIdx, k)) < KAT::magnitude(
zero)) {
781 }
else if (KAT::real(
Rsum(rowIdx, k)) > KAT::magnitude(1.00001)) {
787 for (
size_t k = 0; k < (size_t)
nPDEs; k++) {
793 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) {
794 if (KAT::real(values(entryIdx)) > KAT::real(
zero)) {
795 values(entryIdx) = values(entryIdx) +
801 for (
size_t k = 0; k < (size_t)
nPDEs; k++)
Rsum(rowIdx, k) =
zero;
802 for (
size_t k = 0; k < (size_t)
nPDEs; k++)
nPositive(rowIdx, k) = 0;
807 template <
typename local_matrix_type>
809 using Scalar =
typename local_matrix_type::non_const_value_type;
811 using LO =
typename local_matrix_type::non_const_ordinal_type;
812 using Device =
typename local_matrix_type::device_type;
813 using KAT = Kokkos::ArithTraits<SC>;
820 Kokkos::View<LO**, Device>
inds;
825 origSorted = Kokkos::View<SC**, Device>(
"origSorted", localP_.numRows(), maxRowEntries_);
826 fixedSorted = Kokkos::View<SC**, Device>(
"fixedSorted", localP_.numRows(), maxRowEntries_);
827 inds = Kokkos::View<LO**, Device>(
"inds", localP_.numRows(), maxRowEntries_);
830 KOKKOS_INLINE_FUNCTION
832 auto rowPtr =
localP.graph.row_map;
833 auto values =
localP.values;
835 auto nnz = rowPtr(rowIdx + 1) - rowPtr(rowIdx);
839 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) rsumTarget += values(entryIdx);
843 SC rowSumDeviation, temp, delta;
844 SC closestToLeftBoundDist, closestToRghtBoundDist;
845 LO closestToLeftBound, closestToRghtBound;
850 if ((KAT::real(rsumTarget) >= KAT::real(leftBound * (static_cast<SC>(nnz)))) &&
851 (KAT::real(rsumTarget) <= KAT::real(rghtBound * (static_cast<SC>(nnz))))) {
856 aBigNumber = KAT::zero();
857 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) {
858 if (KAT::magnitude(values(entryIdx)) > KAT::magnitude(aBigNumber))
859 aBigNumber = KAT::magnitude(values(entryIdx));
861 aBigNumber = aBigNumber + (KAT::magnitude(leftBound) + KAT::magnitude(rghtBound)) * (100 *
one);
864 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) {
866 inds(rowIdx, ind) = ind;
870 auto sortVals = Kokkos::subview(
origSorted, rowIdx, Kokkos::ALL());
871 auto sortInds = Kokkos::subview(
inds, rowIdx, Kokkos::make_pair(0, ind));
875 for (
LO i = 1; i < static_cast<LO>(nnz); ++i) {
879 if (KAT::real(sortVals(sortInds(i))) < KAT::real(sortVals(sortInds(0)))) {
880 for (; j > 0; --j) sortInds(j) = sortInds(j - 1);
884 for (; KAT::real(sortVals(ind)) < KAT::real(sortVals(sortInds(j - 1))); --j) sortInds(j) = sortInds(j - 1);
890 for (
LO i = 0; i < static_cast<LO>(nnz); i++)
origSorted(rowIdx, i) = values(rowPtr(rowIdx) +
inds(rowIdx, i));
892 closestToLeftBound = 0;
893 while ((closestToLeftBound < static_cast<LO>(nnz)) && (KAT::real(
origSorted(rowIdx, closestToLeftBound)) <= KAT::real(leftBound))) closestToLeftBound++;
896 closestToRghtBound = closestToLeftBound;
897 while ((closestToRghtBound < static_cast<LO>(nnz)) && (KAT::real(
origSorted(rowIdx, closestToRghtBound)) <= KAT::real(rghtBound))) closestToRghtBound++;
902 closestToLeftBoundDist =
origSorted(rowIdx, closestToLeftBound) - leftBound;
903 if (closestToRghtBound == static_cast<LO>(nnz))
904 closestToRghtBoundDist = aBigNumber;
906 closestToRghtBoundDist =
origSorted(rowIdx, closestToRghtBound) - rghtBound;
911 rowSumDeviation = leftBound * (
static_cast<SC>(closestToLeftBound)) + (
static_cast<SC>(nnz - closestToRghtBound)) * rghtBound - rsumTarget;
912 for (
LO i = closestToLeftBound; i < closestToRghtBound; i++) rowSumDeviation +=
origSorted(rowIdx, i);
917 if (KAT::real(rowSumDeviation) < KAT::real(KAT::zero())) {
920 leftBound = -rghtBound;
927 for (
LO i = 0; i < static_cast<LO>(nnz / 2); i++) {
935 LO itemp = closestToLeftBound;
936 closestToLeftBound = nnz - closestToRghtBound;
937 closestToRghtBound = nnz - itemp;
938 closestToLeftBoundDist =
origSorted(rowIdx, closestToLeftBound) - leftBound;
939 if (closestToRghtBound == static_cast<LO>(nnz))
940 closestToRghtBoundDist = aBigNumber;
942 closestToRghtBoundDist =
origSorted(rowIdx, closestToRghtBound) - rghtBound;
944 rowSumDeviation = -rowSumDeviation;
949 for (
LO i = 0; i < closestToLeftBound; i++)
fixedSorted(rowIdx, i) = leftBound;
951 for (
LO i = closestToRghtBound; i < static_cast<LO>(nnz); i++)
fixedSorted(rowIdx, i) = rghtBound;
953 while ((KAT::magnitude(rowSumDeviation) > KAT::magnitude((
one * 1.e-10) * rsumTarget))) {
954 if (closestToRghtBound != closestToLeftBound)
955 delta = rowSumDeviation /
static_cast<SC>(closestToRghtBound - closestToLeftBound);
959 if (KAT::magnitude(closestToLeftBoundDist) <= KAT::magnitude(closestToRghtBoundDist)) {
960 if (KAT::magnitude(delta) <= KAT::magnitude(closestToLeftBoundDist)) {
961 rowSumDeviation =
zero;
962 for (
LO i = closestToLeftBound; i < closestToRghtBound; i++)
fixedSorted(rowIdx, i) =
origSorted(rowIdx, i) - delta;
964 rowSumDeviation = rowSumDeviation - closestToLeftBoundDist;
965 fixedSorted(rowIdx, closestToLeftBound) = leftBound;
966 closestToLeftBound++;
967 if (closestToLeftBound < static_cast<LO>(nnz))
968 closestToLeftBoundDist =
origSorted(rowIdx, closestToLeftBound) - leftBound;
970 closestToLeftBoundDist = aBigNumber;
973 if (KAT::magnitude(delta) <= KAT::magnitude(closestToRghtBoundDist)) {
975 for (
LO i = closestToLeftBound; i < closestToRghtBound; i++)
fixedSorted(rowIdx, i) =
origSorted(rowIdx, i) - delta;
977 rowSumDeviation = rowSumDeviation + closestToRghtBoundDist;
980 closestToRghtBound++;
982 if (closestToRghtBound >= static_cast<LO>(nnz))
983 closestToRghtBoundDist = aBigNumber;
985 closestToRghtBoundDist =
origSorted(rowIdx, closestToRghtBound) - rghtBound;
990 auto rowStart = rowPtr(rowIdx);
995 for (
LO i = 0; i < static_cast<LO>(nnz / 2); i++) {
1002 for (
LO i = 0; i < static_cast<LO>(nnz); i++) values(rowStart +
inds(rowIdx, i)) =
fixedSorted(rowIdx, i);
1006 for (
auto entryIdx = rowPtr(rowIdx); entryIdx < rowPtr(rowIdx + 1); entryIdx++) values(entryIdx) =
one / (
static_cast<SC>(nnz));
1013 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1015 using Device =
typename Matrix::local_matrix_type::device_type;
1016 LO nPDEs = A->GetFixedBlockSize();
1018 using local_mat_type =
typename Matrix::local_matrix_type;
1020 Kokkos::parallel_for(
"enforce constraint", Kokkos::RangePolicy<typename Device::execution_space>(0, P->getRowMap()->getLocalNumElements()),
1025 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1027 using Device =
typename Matrix::local_matrix_type::device_type;
1030 using local_mat_type =
typename Matrix::local_matrix_type;
1032 Kokkos::parallel_for(
"enforce constraint", Kokkos::RangePolicy<typename Device::execution_space>(0, P->getLocalNumRows()),
1039 #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
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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
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)
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.