119 const std::string prefix =
"MueLu::SaPFactory(" + levelIDs +
"): ";
121 typedef typename Teuchos::ScalarTraits<SC>::coordinateType Coordinate;
122 typedef typename Teuchos::ScalarTraits<SC>::magnitudeType MT;
128 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
129 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.
GetFactoryManager()->GetFactory(
"Ptent"); }
130 const ParameterList& pL = GetParameterList();
133 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
134 RCP<Matrix> Ptent = coarseLevel.
Get< RCP<Matrix> >(
"P", initialPFact.get());
138 if (Ptent == Teuchos::null) {
139 finalP = Teuchos::null;
140 Set(coarseLevel,
"P", finalP);
144 if (restrictionMode_) {
152 RCP<ParameterList> APparams = rcp(
new ParameterList);;
153 if(pL.isSublist(
"matrixmatrix: kernel params"))
154 APparams->sublist(
"matrixmatrix: kernel params") = pL.sublist(
"matrixmatrix: kernel params");
156 if (coarseLevel.
IsAvailable(
"AP reuse data",
this)) {
157 GetOStream(
static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
159 APparams = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data",
this);
161 if (APparams->isParameter(
"graph"))
162 finalP = APparams->get< RCP<Matrix> >(
"graph");
165 APparams->set(
"compute global constants: temporaries",APparams->get(
"compute global constants: temporaries",
false));
166 APparams->set(
"compute global constants",APparams->get(
"compute global constants",
false));
168 const SC
dampingFactor = as<SC>(pL.get<
double>(
"sa: damping factor"));
169 const LO maxEigenIterations = as<LO>(pL.get<
int> (
"sa: eigenvalue estimate num iterations"));
170 const bool estimateMaxEigen = pL.get<
bool> (
"sa: calculate eigenvalue estimate");
171 const bool useAbsValueRowSum= pL.get<
bool> (
"sa: use rowsumabs diagonal scaling");
172 const bool doQRStep = pL.get<
bool>(
"tentative: calculate qr");
173 const bool enforceConstraints= pL.get<
bool>(
"sa: enforce constraints");
174 const MT userDefinedMaxEigen = as<MT>(pL.get<
double>(
"sa: max eigenvalue"));
175 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
176 double dTol = pL.get<
double>(
"sa: rowsumabs diagonal replacement tolerance");
177 const Magnitude diagonalReplacementTolerance = (dTol == as<double>(-1) ? Teuchos::ScalarTraits<Scalar>::eps()*100 : as<Magnitude>(pL.get<
double>(
"sa: rowsumabs diagonal replacement tolerance")));
178 const SC diagonalReplacementValue = as<SC>(pL.get<
double>(
"sa: rowsumabs diagonal replacement value"));
179 const bool replaceSingleEntryRowWithZero = pL.get<
bool>(
"sa: rowsumabs replace single entry row with zero");
180 const bool useAutomaticDiagTol = pL.get<
bool>(
"sa: rowsumabs use automatic diagonal tolerance");
184 "MueLu::TentativePFactory::MakeTentative: cannot use 'enforce constraints' and 'calculate qr' at the same time");
189 Teuchos::RCP<Vector> invDiag;
190 if (userDefinedMaxEigen == -1.)
193 lambdaMax = A->GetMaxEigenvalueEstimate();
194 if (lambdaMax == -Teuchos::ScalarTraits<SC>::one() || estimateMaxEigen) {
195 GetOStream(
Statistics1) <<
"Calculating max eigenvalue estimate now (max iters = "<< maxEigenIterations <<
196 ( (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;
216 lambdaMax = userDefinedMaxEigen;
217 A->SetMaxEigenvalueEstimate(lambdaMax);
223 if (!useAbsValueRowSum)
225 else if (invDiag == Teuchos::null) {
226 GetOStream(
Runtime0) <<
"Using rowsumabs diagonal" << std::endl;
227 const bool returnReciprocal=
true;
229 diagonalReplacementTolerance,
230 diagonalReplacementValue,
231 replaceSingleEntryRowWithZero,
232 useAutomaticDiagTol);
233 TEUCHOS_TEST_FOR_EXCEPTION(invDiag.is_null(),
Exceptions::RuntimeError,
"SaPFactory: diagonal reciprocal is null.");
237 TEUCHOS_TEST_FOR_EXCEPTION(!std::isfinite(Teuchos::ScalarTraits<SC>::magnitude(omega)),
Exceptions::RuntimeError,
"Prolongator damping factor needs to be finite.");
240 finalP = Xpetra::IteratorOps<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Jacobi(omega, *invDiag, *A, *Ptent, finalP,
241 GetOStream(
Statistics2), std::string(
"MueLu::SaP-")+levelIDs, APparams);
242 if (enforceConstraints) {
243 if (A->GetFixedBlockSize() == 1) optimalSatisfyPConstraintsForScalarPDEs( finalP);
244 else SatisfyPConstraints( A, finalP);
253 if (!restrictionMode_) {
255 if(!finalP.is_null()) {std::ostringstream oss; oss <<
"P_" << coarseLevel.
GetLevelID(); finalP->setObjectLabel(oss.str());}
256 Set(coarseLevel,
"P", finalP);
258 APparams->set(
"graph", finalP);
259 Set(coarseLevel,
"AP reuse data", APparams);
262 if (Ptent->IsView(
"stridedMaps"))
263 finalP->CreateView(
"stridedMaps", Ptent);
266 RCP<ParameterList> params = rcp(
new ParameterList());
267 params->set(
"printLoadBalancingInfo",
true);
268 params->set(
"printCommInfo",
true);
278 if(!R.is_null()) {std::ostringstream oss; oss <<
"R_" << coarseLevel.
GetLevelID(); R->setObjectLabel(oss.str());}
281 Set(coarseLevel,
"R", R);
284 if (Ptent->IsView(
"stridedMaps"))
285 R->CreateView(
"stridedMaps", Ptent,
true);
288 RCP<ParameterList> params = rcp(
new ParameterList());
289 params->set(
"printLoadBalancingInfo",
true);
290 params->set(
"printCommInfo",
true);
316 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
317 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
318 LO nPDEs = A->GetFixedBlockSize();
319 Teuchos::ArrayRCP<Scalar> ConstraintViolationSum(nPDEs);
320 Teuchos::ArrayRCP<Scalar> Rsum(nPDEs);
321 Teuchos::ArrayRCP<size_t> nPositive(nPDEs);
322 for (
size_t k=0; k < (size_t) nPDEs; k++) ConstraintViolationSum[k] = zero;
323 for (
size_t k=0; k < (size_t) nPDEs; k++) Rsum[k] = zero;
324 for (
size_t k=0; k < (size_t) nPDEs; k++) nPositive[k] = 0;
327 for (
size_t i = 0; i < as<size_t>(P->getRowMap()->getLocalNumElements()); i++) {
329 Teuchos::ArrayView<const LocalOrdinal> indices;
330 Teuchos::ArrayView<const Scalar> vals1;
331 Teuchos::ArrayView< Scalar> vals;
332 P->getLocalRowView((LO) i, indices, vals1);
333 size_t nnz = indices.size();
334 if (nnz == 0)
continue;
336 vals = ArrayView<Scalar>(
const_cast<SC*
>(vals1.getRawPtr()), nnz);
339 bool checkRow =
true;
345 for (LO j = 0; j < indices.size(); j++) {
346 Rsum[ j%nPDEs ] += vals[j];
347 if (Teuchos::ScalarTraits<SC>::real(vals[j]) < Teuchos::ScalarTraits<SC>::real(zero)) {
348 ConstraintViolationSum[ j%nPDEs ] += vals[j];
352 if (Teuchos::ScalarTraits<SC>::real(vals[j]) != Teuchos::ScalarTraits<SC>::real(zero))
353 (nPositive[ j%nPDEs])++;
355 if (Teuchos::ScalarTraits<SC>::real(vals[j]) > Teuchos::ScalarTraits<SC>::real(1.00001 )) {
356 ConstraintViolationSum[ j%nPDEs ] += (vals[j] - one);
366 for (
size_t k=0; k < (size_t) nPDEs; k++) {
368 if (Teuchos::ScalarTraits<SC>::real(Rsum[ k ]) < Teuchos::ScalarTraits<SC>::magnitude(zero)) {
369 ConstraintViolationSum[k] += (-Rsum[k]);
371 else if (Teuchos::ScalarTraits<SC>::real(Rsum[ k ]) > Teuchos::ScalarTraits<SC>::magnitude(1.00001)) {
372 ConstraintViolationSum[k] += (one - Rsum[k]);
377 for (
size_t k=0; k < (size_t) nPDEs; k++) {
378 if (Teuchos::ScalarTraits<SC>::magnitude(ConstraintViolationSum[ k ]) != Teuchos::ScalarTraits<SC>::magnitude(zero))
383 for (LO j = 0; j < indices.size(); j++) {
384 if (Teuchos::ScalarTraits<SC>::real(vals[j]) > Teuchos::ScalarTraits<SC>::real(zero)) {
385 vals[j] += (ConstraintViolationSum[j%nPDEs]/ (as<Scalar>(nPositive[j%nPDEs])));
388 for (
size_t k=0; k < (size_t) nPDEs; k++) ConstraintViolationSum[k] = zero;
390 for (
size_t k=0; k < (size_t) nPDEs; k++) Rsum[k] = zero;
391 for (
size_t k=0; k < (size_t) nPDEs; k++) nPositive[k] = 0;
399 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
400 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
403 Teuchos::ArrayRCP<Scalar> scalarData(3*maxEntriesPerRow);
406 for (
size_t i = 0; i < as<size_t>(P->getRowMap()->getLocalNumElements()); i++) {
408 Teuchos::ArrayView<const LocalOrdinal> indices;
409 Teuchos::ArrayView<const Scalar> vals1;
410 Teuchos::ArrayView< Scalar> vals;
412 size_t nnz = indices.size();
415 vals = ArrayView<Scalar>(
const_cast<SC*
>(vals1.getRawPtr()), nnz);
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);
552 Scalar notFlippedLeftBound, notFlippedRghtBound, aBigNumber, *origSorted;
553 Scalar rowSumDeviation, temp, *fixedSorted, delta;
554 Scalar closestToLeftBoundDist, closestToRghtBoundDist;
559 notFlippedLeftBound = leftBound;
560 notFlippedRghtBound = rghtBound;
562 if ((Teuchos::ScalarTraits<SC>::real(rsumTarget) >= Teuchos::ScalarTraits<SC>::real(leftBound*as<Scalar>(nEntries))) &&
563 (Teuchos::ScalarTraits<SC>::real(rsumTarget) <= Teuchos::ScalarTraits<SC>::real(rghtBound*as<Scalar>(nEntries))))
564 hasFeasibleSol =
true;
566 hasFeasibleSol=
false;
567 return hasFeasibleSol;
572 aBigNumber = Teuchos::ScalarTraits<SC>::zero();
574 if ( Teuchos::ScalarTraits<SC>::magnitude(orig[i]) > Teuchos::ScalarTraits<SC>::magnitude(aBigNumber))
575 aBigNumber = Teuchos::ScalarTraits<SC>::magnitude(orig[i]);
577 aBigNumber = aBigNumber+ (Teuchos::ScalarTraits<SC>::magnitude(leftBound) + Teuchos::ScalarTraits<SC>::magnitude(rghtBound))*as<Scalar>(100.0);
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,
589 {
return Teuchos::ScalarTraits<SC>::real(origSorted[leftIndex]) < Teuchos::ScalarTraits<SC>::real(origSorted[rightIndex]);});
591 for (
LocalOrdinal i = 0; i < nEntries; i++) origSorted[i] = orig[inds[i]];
593 closestToLeftBound = 0;
594 while ((closestToLeftBound < nEntries) && (Teuchos::ScalarTraits<SC>::real(origSorted[closestToLeftBound]) <= Teuchos::ScalarTraits<SC>::real(leftBound))) closestToLeftBound++;
597 closestToRghtBound = closestToLeftBound;
598 while ((closestToRghtBound < nEntries) && (Teuchos::ScalarTraits<SC>::real(origSorted[closestToRghtBound]) <= Teuchos::ScalarTraits<SC>::real(rghtBound))) closestToRghtBound++;
603 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
604 if (closestToRghtBound==nEntries) closestToRghtBoundDist= aBigNumber;
605 else closestToRghtBoundDist= origSorted[closestToRghtBound] - rghtBound;
610 rowSumDeviation = leftBound*as<Scalar>(closestToLeftBound) + as<Scalar>((nEntries-closestToRghtBound))*rghtBound - rsumTarget;
611 for (
LocalOrdinal i=closestToLeftBound; i < closestToRghtBound; i++) rowSumDeviation += origSorted[i];
616 if (Teuchos::ScalarTraits<SC>::real(rowSumDeviation) < Teuchos::ScalarTraits<SC>::real(Teuchos::ScalarTraits<SC>::zero())) {
618 temp = leftBound; leftBound = -rghtBound; rghtBound = temp;
622 if ((nEntries%2) == 1) origSorted[(nEntries/2) ] = -origSorted[(nEntries/2) ];
625 origSorted[i] = -origSorted[nEntries-1-i];
626 origSorted[nEntries-i-1] = -temp;
632 closestToLeftBound = nEntries-closestToRghtBound;
633 closestToRghtBound = nEntries-itemp;
634 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
635 if (closestToRghtBound==nEntries) closestToRghtBoundDist= aBigNumber;
636 else closestToRghtBoundDist= origSorted[closestToRghtBound] - rghtBound;
638 rowSumDeviation = -rowSumDeviation;
643 for (
LocalOrdinal i = 0; i < closestToLeftBound; i++) fixedSorted[i] = leftBound;
644 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound; i++) fixedSorted[i] = origSorted[i];
645 for (
LocalOrdinal i = closestToRghtBound; i < nEntries; i++) fixedSorted[i] = rghtBound;
647 while ((Teuchos::ScalarTraits<SC>::magnitude(rowSumDeviation) > Teuchos::ScalarTraits<SC>::magnitude(as<Scalar>(1.e-10)*rsumTarget))){
648 if (closestToRghtBound != closestToLeftBound)
649 delta = rowSumDeviation/ as<Scalar>(closestToRghtBound - closestToLeftBound);
650 else delta = aBigNumber;
652 if (Teuchos::ScalarTraits<SC>::magnitude(closestToLeftBoundDist) <= Teuchos::ScalarTraits<SC>::magnitude(closestToRghtBoundDist)) {
653 if (Teuchos::ScalarTraits<SC>::magnitude(delta) <= Teuchos::ScalarTraits<SC>::magnitude(closestToLeftBoundDist)) {
654 rowSumDeviation = Teuchos::ScalarTraits<SC>::zero();
655 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound ; i++) fixedSorted[i] = origSorted[i] - delta;
658 rowSumDeviation = rowSumDeviation - closestToLeftBoundDist;
659 fixedSorted[closestToLeftBound] = leftBound;
660 closestToLeftBound++;
661 if (closestToLeftBound < nEntries) closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
662 else closestToLeftBoundDist = aBigNumber;
666 if (Teuchos::ScalarTraits<SC>::magnitude(delta) <= Teuchos::ScalarTraits<SC>::magnitude(closestToRghtBoundDist)) {
668 for (
LocalOrdinal i = closestToLeftBound; i < closestToRghtBound ; i++) fixedSorted[i] = origSorted[i] - delta;
671 rowSumDeviation = rowSumDeviation + closestToRghtBoundDist;
673 fixedSorted[closestToRghtBound] = origSorted[closestToRghtBound];
674 closestToRghtBound++;
676 if (closestToRghtBound >= nEntries) closestToRghtBoundDist = aBigNumber;
677 else closestToRghtBoundDist = origSorted[closestToRghtBound] - rghtBound;
685 if ((nEntries%2) == 1) fixedSorted[(nEntries/2) ] = -fixedSorted[(nEntries/2) ];
688 fixedSorted[i] = -fixedSorted[nEntries-1-i];
689 fixedSorted[nEntries-i-1] = -temp;
692 for (
LocalOrdinal i = 0; i < nEntries; i++) fixedUnsorted[inds[i]] = fixedSorted[i];
696 bool lowerViolation =
false;
697 bool upperViolation =
false;
698 bool sumViolation =
false;
699 using TST = Teuchos::ScalarTraits<SC>;
702 if (TST::real(fixedUnsorted[i]) < TST::real(notFlippedLeftBound)) lowerViolation =
true;
703 if (TST::real(fixedUnsorted[i]) > TST::real(notFlippedRghtBound)) upperViolation =
true;
704 temp += fixedUnsorted[i];
706 SC tol = as<Scalar>(std::max(1.0e-8, as<double>(100*TST::eps())));
707 if (TST::magnitude(temp - rsumTarget) > TST::magnitude(tol*rsumTarget)) sumViolation =
true;
709 TEUCHOS_TEST_FOR_EXCEPTION(lowerViolation,
Exceptions::RuntimeError,
"MueLu::SaPFactory::constrainRow: feasible solution but computation resulted in a lower bound violation??? ");
710 TEUCHOS_TEST_FOR_EXCEPTION(upperViolation,
Exceptions::RuntimeError,
"MueLu::SaPFactory::constrainRow: feasible solution but computation resulted in an upper bound violation??? ");
711 TEUCHOS_TEST_FOR_EXCEPTION(sumViolation,
Exceptions::RuntimeError,
"MueLu::SaPFactory::constrainRow: feasible solution but computation resulted in a row sum violation??? ");
713 return hasFeasibleSol;