MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_SaPFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_SAPFACTORY_DEF_HPP
47#define MUELU_SAPFACTORY_DEF_HPP
48
49#include <Xpetra_Matrix.hpp>
50#include <Xpetra_IteratorOps.hpp> // containing routines to generate Jacobi iterator
51#include <Xpetra_IO.hpp>
52#include <sstream>
53
55
57#include "MueLu_Level.hpp"
58#include "MueLu_MasterList.hpp"
59#include "MueLu_Monitor.hpp"
60#include "MueLu_PerfUtils.hpp"
62#include "MueLu_TentativePFactory.hpp"
63#include "MueLu_Utilities.hpp"
64
65namespace MueLu {
66
67 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69 RCP<ParameterList> validParamList = rcp(new ParameterList());
70
71#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
72 SET_VALID_ENTRY("sa: damping factor");
73 SET_VALID_ENTRY("sa: calculate eigenvalue estimate");
74 SET_VALID_ENTRY("sa: eigenvalue estimate num iterations");
75 SET_VALID_ENTRY("sa: use rowsumabs diagonal scaling");
76 SET_VALID_ENTRY("sa: enforce constraints");
77 SET_VALID_ENTRY("tentative: calculate qr");
78 SET_VALID_ENTRY("sa: max eigenvalue");
79 SET_VALID_ENTRY("sa: rowsumabs diagonal replacement tolerance");
80 SET_VALID_ENTRY("sa: rowsumabs diagonal replacement value");
81 SET_VALID_ENTRY("sa: rowsumabs replace single entry row with zero");
82 SET_VALID_ENTRY("sa: rowsumabs use automatic diagonal tolerance");
83#undef SET_VALID_ENTRY
84
85 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A used during the prolongator smoothing process");
86 validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Tentative prolongator factory");
87
88 // Make sure we don't recursively validate options for the matrixmatrix kernels
89 ParameterList norecurse;
90 norecurse.disableRecursiveValidation();
91 validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
92
93
94 return validParamList;
95 }
96
97 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
99 Input(fineLevel, "A");
100
101 // Get default tentative prolongator factory
102 // Getting it that way ensure that the same factory instance will be used for both SaPFactory and NullspaceFactory.
103 RCP<const FactoryBase> initialPFact = GetFactory("P");
104 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent"); }
105 coarseLevel.DeclareInput("P", initialPFact.get(), this); // --
106 }
107
108 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
110 return BuildP(fineLevel, coarseLevel);
111 }
112
113 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
115 FactoryMonitor m(*this, "Prolongator smoothing", coarseLevel);
116
117 std::string levelIDs = toString(coarseLevel.GetLevelID());
118
119 const std::string prefix = "MueLu::SaPFactory(" + levelIDs + "): ";
120
121 typedef typename Teuchos::ScalarTraits<SC>::coordinateType Coordinate;
122 typedef typename Teuchos::ScalarTraits<SC>::magnitudeType MT;
123
124
125 // Get default tentative prolongator factory
126 // Getting it that way ensure that the same factory instance will be used for both SaPFactory and NullspaceFactory.
127 // -- Warning: Do not use directly initialPFact_. Use initialPFact instead everywhere!
128 RCP<const FactoryBase> initialPFact = GetFactory("P");
129 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.GetFactoryManager()->GetFactory("Ptent"); }
130 const ParameterList& pL = GetParameterList();
131
132 // Level Get
133 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
134 RCP<Matrix> Ptent = coarseLevel.Get< RCP<Matrix> >("P", initialPFact.get());
135 RCP<Matrix> finalP;
136 // If Tentative facctory bailed out (e.g., number of global aggregates is 0), then SaPFactory bails
137 // This level will ultimately be removed in MueLu_Hierarchy_defs.h via a resize()
138 if (Ptent == Teuchos::null) {
139 finalP = Teuchos::null;
140 Set(coarseLevel, "P", finalP);
141 return;
142 }
143
144 if (restrictionMode_) {
145 SubFactoryMonitor m2(*this, "Transpose A", coarseLevel);
146 A = Utilities::Transpose(*A, true); // build transpose of A explicitely
147 }
148
149 // Build final prolongator
150
151 // Reuse pattern if available
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");
155
156 if (coarseLevel.IsAvailable("AP reuse data", this)) {
157 GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous AP data" << std::endl;
158
159 APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
160
161 if (APparams->isParameter("graph"))
162 finalP = APparams->get< RCP<Matrix> >("graph");
163 }
164 // By default, we don't need global constants for SaP
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));
167
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");
181
182 // Sanity checking
183 TEUCHOS_TEST_FOR_EXCEPTION(doQRStep && enforceConstraints,Exceptions::RuntimeError,
184 "MueLu::TentativePFactory::MakeTentative: cannot use 'enforce constraints' and 'calculate qr' at the same time");
185
186 if (dampingFactor != Teuchos::ScalarTraits<SC>::zero()) {
187
188 Scalar lambdaMax;
189 Teuchos::RCP<Vector> invDiag;
190 if (userDefinedMaxEigen == -1.)
191 {
192 SubFactoryMonitor m2(*this, "Eigenvalue estimate", coarseLevel);
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;
200 invDiag = Utilities::GetLumpedMatrixDiagonal(*A,returnReciprocal,
201 diagonalReplacementTolerance,
202 diagonalReplacementValue,
203 replaceSingleEntryRowWithZero,
204 useAutomaticDiagTol);
205 TEUCHOS_TEST_FOR_EXCEPTION(invDiag.is_null(), Exceptions::RuntimeError,
206 "SaPFactory: eigenvalue estimate: diagonal reciprocal is null.");
207 lambdaMax = Utilities::PowerMethod(*A, invDiag, maxEigenIterations, stopTol);
208 } else
209 lambdaMax = Utilities::PowerMethod(*A, true, maxEigenIterations, stopTol);
210 A->SetMaxEigenvalueEstimate(lambdaMax);
211 } else {
212 GetOStream(Statistics1) << "Using cached max eigenvalue estimate" << std::endl;
213 }
214 }
215 else {
216 lambdaMax = userDefinedMaxEigen;
217 A->SetMaxEigenvalueEstimate(lambdaMax);
218 }
219 GetOStream(Statistics1) << "Prolongator damping factor = " << dampingFactor/lambdaMax << " (" << dampingFactor << " / " << lambdaMax << ")" << std::endl;
220
221 {
222 SubFactoryMonitor m2(*this, "Fused (I-omega*D^{-1} A)*Ptent", coarseLevel);
223 if (!useAbsValueRowSum)
224 invDiag = Utilities::GetMatrixDiagonalInverse(*A); //default
225 else if (invDiag == Teuchos::null) {
226 GetOStream(Runtime0) << "Using rowsumabs diagonal" << std::endl;
227 const bool returnReciprocal=true;
228 invDiag = Utilities::GetLumpedMatrixDiagonal(*A,returnReciprocal,
229 diagonalReplacementTolerance,
230 diagonalReplacementValue,
231 replaceSingleEntryRowWithZero,
232 useAutomaticDiagTol);
233 TEUCHOS_TEST_FOR_EXCEPTION(invDiag.is_null(), Exceptions::RuntimeError, "SaPFactory: diagonal reciprocal is null.");
234 }
235
236 SC omega = dampingFactor / lambdaMax;
237 TEUCHOS_TEST_FOR_EXCEPTION(!std::isfinite(Teuchos::ScalarTraits<SC>::magnitude(omega)), Exceptions::RuntimeError, "Prolongator damping factor needs to be finite.");
238
239 // finalP = Ptent + (I - \omega D^{-1}A) Ptent
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);
245 }
246 }
247
248 } else {
249 finalP = Ptent;
250 }
251
252 // Level Set
253 if (!restrictionMode_) {
254 // The factory is in prolongation mode
255 if(!finalP.is_null()) {std::ostringstream oss; oss << "P_" << coarseLevel.GetLevelID(); finalP->setObjectLabel(oss.str());}
256 Set(coarseLevel, "P", finalP);
257
258 APparams->set("graph", finalP);
259 Set(coarseLevel, "AP reuse data", APparams);
260
261 // NOTE: EXPERIMENTAL
262 if (Ptent->IsView("stridedMaps"))
263 finalP->CreateView("stridedMaps", Ptent);
264
265 if (IsPrint(Statistics2)) {
266 RCP<ParameterList> params = rcp(new ParameterList());
267 params->set("printLoadBalancingInfo", true);
268 params->set("printCommInfo", true);
269 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*finalP, "P", params);
270 }
271
272 } else {
273 // The factory is in restriction mode
274 RCP<Matrix> R;
275 {
276 SubFactoryMonitor m2(*this, "Transpose P", coarseLevel);
277 R = Utilities::Transpose(*finalP, true);
278 if(!R.is_null()) {std::ostringstream oss; oss << "R_" << coarseLevel.GetLevelID(); R->setObjectLabel(oss.str());}
279 }
280
281 Set(coarseLevel, "R", R);
282
283 // NOTE: EXPERIMENTAL
284 if (Ptent->IsView("stridedMaps"))
285 R->CreateView("stridedMaps", Ptent, true/*transposeA*/);
286
287 if (IsPrint(Statistics2)) {
288 RCP<ParameterList> params = rcp(new ParameterList());
289 params->set("printLoadBalancingInfo", true);
290 params->set("printCommInfo", true);
291 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*R, "R", params);
292 }
293 }
294
295 } //Build()
296
297 // Analyze the grid transfer produced by smoothed aggregation and make
298 // modifications if it does not look right. In particular, if there are
299 // negative entries or entries larger than 1, modify P's rows.
300 //
301 // Note: this kind of evaluation probably only makes sense if not doing QR
302 // when constructing tentative P.
303 //
304 // For entries that do not satisfy the two constraints (>= 0 or <=1) we set
305 // these entries to the constraint value and modify the rest of the row
306 // so that the row sum remains the same as before by adding an equal
307 // amount to each remaining entry. However, if the original row sum value
308 // violates the constraints, we set the row sum back to 1 (the row sum of
309 // tentative P). After doing the modification to a row, we need to check
310 // again the entire row to make sure that the modified row does not violate
311 // the constraints.
312
313 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
315
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;
325
326
327 for (size_t i = 0; i < as<size_t>(P->getRowMap()->getLocalNumElements()); i++) {
328
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;
335
336 vals = ArrayView<Scalar>(const_cast<SC*>(vals1.getRawPtr()), nnz);
337
338
339 bool checkRow = true;
340
341 while (checkRow) {
342
343 // check constraints and compute the row sum
344
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];
349 vals[j] = zero;
350 }
351 else {
352 if (Teuchos::ScalarTraits<SC>::real(vals[j]) != Teuchos::ScalarTraits<SC>::real(zero))
353 (nPositive[ j%nPDEs])++;
354
355 if (Teuchos::ScalarTraits<SC>::real(vals[j]) > Teuchos::ScalarTraits<SC>::real(1.00001 )) {
356 ConstraintViolationSum[ j%nPDEs ] += (vals[j] - one);
357 vals[j] = one;
358 }
359 }
360 }
361
362 checkRow = false;
363
364 // take into account any row sum that violates the contraints
365
366 for (size_t k=0; k < (size_t) nPDEs; k++) {
367
368 if (Teuchos::ScalarTraits<SC>::real(Rsum[ k ]) < Teuchos::ScalarTraits<SC>::magnitude(zero)) {
369 ConstraintViolationSum[k] += (-Rsum[k]); // rstumin
370 }
371 else if (Teuchos::ScalarTraits<SC>::real(Rsum[ k ]) > Teuchos::ScalarTraits<SC>::magnitude(1.00001)) {
372 ConstraintViolationSum[k] += (one - Rsum[k]); // rstumin
373 }
374 }
375
376 // check if row need modification
377 for (size_t k=0; k < (size_t) nPDEs; k++) {
378 if (Teuchos::ScalarTraits<SC>::magnitude(ConstraintViolationSum[ k ]) != Teuchos::ScalarTraits<SC>::magnitude(zero))
379 checkRow = true;
380 }
381 // modify row
382 if (checkRow) {
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])));
386 }
387 }
388 for (size_t k=0; k < (size_t) nPDEs; k++) ConstraintViolationSum[k] = zero;
389 }
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;
392 } // while (checkRow) ...
393 } // for (size_t i = 0; i < as<size_t>(P->getRowMap()->getNumNodeElements()); i++) ...
394 } //SatsifyPConstraints()
395
396 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
398
399 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
400 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
401
402 LocalOrdinal maxEntriesPerRow = 100; // increased later if needed
403 Teuchos::ArrayRCP<Scalar> scalarData(3*maxEntriesPerRow);
404 bool hasFeasible;
405
406 for (size_t i = 0; i < as<size_t>(P->getRowMap()->getLocalNumElements()); i++) {
407
408 Teuchos::ArrayView<const LocalOrdinal> indices;
409 Teuchos::ArrayView<const Scalar> vals1;
410 Teuchos::ArrayView< Scalar> vals;
411 P->getLocalRowView((LocalOrdinal) i, indices, vals1);
412 size_t nnz = indices.size();
413 if (nnz != 0) {
414
415 vals = ArrayView<Scalar>(const_cast<SC*>(vals1.getRawPtr()), nnz);
416 Scalar rsumTarget = zero;
417 for (size_t j = 0; j < nnz; j++) rsumTarget += vals[j];
418
419 if (nnz > as<size_t>(maxEntriesPerRow)) {
420 maxEntriesPerRow = nnz*3;
421 scalarData.resize(3*maxEntriesPerRow);
422 }
423 hasFeasible = constrainRow(vals.getRawPtr(), as<LocalOrdinal>(nnz), zero , one, rsumTarget, vals.getRawPtr(), scalarData.getRawPtr());
424
425 if (!hasFeasible) { // just set all entries to the same value giving a row sum of 1
426 for (size_t j = 0; j < nnz; j++) vals[j] = one/as<Scalar>(nnz);
427 }
428 }
429
430 } // for (size_t i = 0; i < as<size_t>(P->getRowMap()->getNumNodeElements()); i++) ...
431 } //SatsifyPConstraints()
432
433 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
434 bool SaPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::constrainRow(Scalar *orig, LocalOrdinal nEntries, Scalar leftBound, Scalar rghtBound,Scalar rsumTarget, Scalar *fixedUnsorted, Scalar *scalarData) const {
435/*
436 Input
437 orig data that should be adjusted to satisfy bound constraints and
438 row sum constraint. orig is not modified by this function.
439
440 nEntries length or 'orig'
441
442 leftBound, define bound constraints for the results.
443 rghtBound
444
445 rsumTarget defines an equality constraint for the row sum
446
447 fixedUnsorted on output, if a feasible solutuion exists then
448 || orig - fixedUnsorted || = min when also
449 leftBound <= fixedUnsorted[i] <= rghtBound for all i
450 and sum(fixedUnsorted) = rsumTarget.
451
452 Note: it is possible to use the same pointer for
453 fixedUnsorted and orig. In this case, orig gets
454 overwritten with the new constraint satisfying values.
455
456 scalarData a work array that should be 3x nEntries.
457
458 On return constrain() indicates whether or not a feasible solution exists.
459*/
460
461/*
462 Given a sequence of numbers o1 ... on, fix these so that they are as
463 close as possible to the original but satisfy bound constraints and also
464 have the same row sum as the oi's. If we know who is going to lie on a
465 bound, then the "best" answer (i.e., || o - f ||_2 = min) perturbs
466 each element that doesn't lie on a bound by the same amount.
467
468 We can represent the oi's by considering scattered points on a number line
469
470 | |
471 | |
472 o o o | o o o o o o |o o
473 | |
474
475 \____/ \____/
476 <---- <----
477 delta delta
478
479 Bounds are shown by vertical lines. The fi's must lie within the bounds, so
480 the 3 leftmost points must be shifted to the right and the 2 rightmost must
481 be shifted to the left. If these fi points are all shifted to the bounds
482 while the others remain in the same location, the row sum constraint is
483 likely not satisfied and so more shifting is necessary. In the figure, the f
484 rowsum is too large and so there must be more shifting to the left.
485
486 To minimize || o - f ||_2, we basically shift all "interiors" by the same
487 amount, denoted delta. The only trick is that some points near bounds are
488 still affected by the bounds and so these points might be shifted more or less
489 than delta. In the example,t he 3 rightmost points are shifted in the opposite
490 direction as delta to the bound. The 4th point is shifted by something less
491 than delta so it does not violate the lower bound. The rightmost point is
492 shifted to the bound by some amount larger than delta. However, the 2nd point
493 is shifted by delta (i.e., it lies inside the two bounds).
494
495 If we know delta, we can figure out everything. If we know which points
496 are special (not shifted by delta), we can also figure out everything.
497 The problem is these two things (delta and the special points) are
498 inter-connected. An algorithm for computing follows.
499
500 1) move exterior points to the bounds and compute how much the row sum is off
501 (rowSumDeviation). We assume now that the new row sum is high, so interior
502 points must be shifted left.
503
504 2) Mark closest point just left of the leftmost bound, closestToLeftBound,
505 and compute its distance to the leftmost bound. Mark closest point to the
506 left of the rightmost bound, closestToRghtBound, and compute its distance to
507 right bound. There are two cases to consider.
508
509 3) Case 1: closestToLeftBound is closer than closestToRghtBound.
510 Assume that shifting by delta does not move closestToLeftBound past the
511 left bound. This means that it will be shifted by delta. However,
512 closestToRghtBound will be shifted by more than delta. So the total
513 number of points shifted by delta (|interiorToBounds|) includes
514 closestToLeftBound up to and including the point just to the left of
515 closestToRghtBound. So
516
517 delta = rowSumDeviation/ |interiorToBounds| .
518
519 Recall that rowSumDeviation already accounts for the non-delta shift of
520 of closestToRightBound. Now check whether our assumption is valid.
521
522 If delta <= closestToLeftBoundDist, assumption is true so delta can be
523 applied to interiorToBounds ... and we are done.
524 Else assumption is false. Shift closestToLeftBound to the left bound.
525 Update rowSumDeviation, interiorToBounds, and identify new
526 closestToLeftBound. Repeat step 3).
527
528 Case 2: closestToRghtBound is closer than closestToLeftBound.
529 Assume that shifting by delta does not move closestToRghtBound past right
530 bound. This means that it must be shifted by more than delta to right
531 bound. So the total number of points shifted by delta again includes
532 closestToLeftBound up to and including the point just to the left of
533 closestToRghtBound. So again compute
534
535 delta = rowSumDeviation/ |interiorToBounds| .
536
537 If delta <= closestToRghtBoundDist, assumption is true so delta is
538 can be applied to interiorToBounds ... and we are done
539 Else assumption is false. Put closestToRghtBound in the
540 interiorToBounds set. Remove it's contribution to rowSumDeviation,
541 identify new closestToRghtBound. Repeat step 3)
542
543
544 To implement, sort the oi's so things like closestToLeftBound is just index
545 into sorted array. Updaing closestToLeftBound or closestToRghtBound requires
546 increment by 1. |interiorToBounds|= closestToRghtBound - closestToLeftBound
547 To handle the case when the rowsum is low (requiring right interior shifts),
548 just flip the signs on data and use the left-shift code (and then flip back
549 before exiting function.
550*/
551 bool hasFeasibleSol;
552 Scalar notFlippedLeftBound, notFlippedRghtBound, aBigNumber, *origSorted;
553 Scalar rowSumDeviation, temp, *fixedSorted, delta;
554 Scalar closestToLeftBoundDist, closestToRghtBoundDist;
555 LocalOrdinal closestToLeftBound, closestToRghtBound;
556 LocalOrdinal *inds;
557 bool flipped;
558
559 notFlippedLeftBound = leftBound;
560 notFlippedRghtBound = rghtBound;
561
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;
565 else {
566 hasFeasibleSol=false;
567 return hasFeasibleSol;
568 }
569 flipped = false;
570 // compute aBigNumber to handle some corner cases where we need
571 // something large so that an if statement will be false
572 aBigNumber = Teuchos::ScalarTraits<SC>::zero();
573 for (LocalOrdinal i = 0; i < nEntries; i++) {
574 if ( Teuchos::ScalarTraits<SC>::magnitude(orig[i]) > Teuchos::ScalarTraits<SC>::magnitude(aBigNumber))
575 aBigNumber = Teuchos::ScalarTraits<SC>::magnitude(orig[i]);
576 }
577 aBigNumber = aBigNumber+ (Teuchos::ScalarTraits<SC>::magnitude(leftBound) + Teuchos::ScalarTraits<SC>::magnitude(rghtBound))*as<Scalar>(100.0);
578
579 origSorted = &scalarData[0];
580 fixedSorted = &(scalarData[nEntries]);
581 inds = (LocalOrdinal *) &(scalarData[2*nEntries]);
582
583 for (LocalOrdinal i = 0; i < nEntries; i++) inds[i] = i;
584 for (LocalOrdinal i = 0; i < nEntries; i++) origSorted[i] = orig[i]; /* orig no longer used */
585
586 // sort so that orig[inds] is sorted.
587 std::sort(inds, inds+nEntries,
588 [origSorted](LocalOrdinal leftIndex, LocalOrdinal rightIndex)
589 { return Teuchos::ScalarTraits<SC>::real(origSorted[leftIndex]) < Teuchos::ScalarTraits<SC>::real(origSorted[rightIndex]);});
590
591 for (LocalOrdinal i = 0; i < nEntries; i++) origSorted[i] = orig[inds[i]];
592 // find entry in origSorted just to the right of the leftBound
593 closestToLeftBound = 0;
594 while ((closestToLeftBound < nEntries) && (Teuchos::ScalarTraits<SC>::real(origSorted[closestToLeftBound]) <= Teuchos::ScalarTraits<SC>::real(leftBound))) closestToLeftBound++;
595
596 // find entry in origSorted just to the right of the rghtBound
597 closestToRghtBound = closestToLeftBound;
598 while ((closestToRghtBound < nEntries) && (Teuchos::ScalarTraits<SC>::real(origSorted[closestToRghtBound]) <= Teuchos::ScalarTraits<SC>::real(rghtBound))) closestToRghtBound++;
599
600 // compute distance between closestToLeftBound and the left bound and the
601 // distance between closestToRghtBound and the right bound.
602
603 closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
604 if (closestToRghtBound==nEntries) closestToRghtBoundDist= aBigNumber;
605 else closestToRghtBoundDist= origSorted[closestToRghtBound] - rghtBound;
606
607 // compute how far the rowSum is off from the target row sum taking into account
608 // numbers that have been shifted to satisfy bound constraint
609
610 rowSumDeviation = leftBound*as<Scalar>(closestToLeftBound) + as<Scalar>((nEntries-closestToRghtBound))*rghtBound - rsumTarget;
611 for (LocalOrdinal i=closestToLeftBound; i < closestToRghtBound; i++) rowSumDeviation += origSorted[i];
612
613 // the code that follow after this if statement assumes that rowSumDeviation is positive. If this
614 // is not the case, flip the signs of everything so that rowSumDeviation is now positive.
615 // Later we will flip the data back to its original form.
616 if (Teuchos::ScalarTraits<SC>::real(rowSumDeviation) < Teuchos::ScalarTraits<SC>::real(Teuchos::ScalarTraits<SC>::zero())) {
617 flipped = true;
618 temp = leftBound; leftBound = -rghtBound; rghtBound = temp;
619
620 /* flip sign of origSorted and reverse ordering so that the negative version is sorted */
621
622 if ((nEntries%2) == 1) origSorted[(nEntries/2) ] = -origSorted[(nEntries/2) ];
623 for (LocalOrdinal i=0; i < nEntries/2; i++) {
624 temp=origSorted[i];
625 origSorted[i] = -origSorted[nEntries-1-i];
626 origSorted[nEntries-i-1] = -temp;
627 }
628
629 /* reverse bounds */
630
631 LocalOrdinal itemp = closestToLeftBound;
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;
637
638 rowSumDeviation = -rowSumDeviation;
639 }
640
641 // initial fixedSorted so bounds are satisfied and interiors correspond to origSorted
642
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;
646
647 while ((Teuchos::ScalarTraits<SC>::magnitude(rowSumDeviation) > Teuchos::ScalarTraits<SC>::magnitude(as<Scalar>(1.e-10)*rsumTarget))){ // && ( (closestToLeftBound < nEntries ) || (closestToRghtBound < nEntries))) {
648 if (closestToRghtBound != closestToLeftBound)
649 delta = rowSumDeviation/ as<Scalar>(closestToRghtBound - closestToLeftBound);
650 else delta = aBigNumber;
651
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;
656 }
657 else {
658 rowSumDeviation = rowSumDeviation - closestToLeftBoundDist;
659 fixedSorted[closestToLeftBound] = leftBound;
660 closestToLeftBound++;
661 if (closestToLeftBound < nEntries) closestToLeftBoundDist = origSorted[closestToLeftBound] - leftBound;
662 else closestToLeftBoundDist = aBigNumber;
663 }
664 }
665 else {
666 if (Teuchos::ScalarTraits<SC>::magnitude(delta) <= Teuchos::ScalarTraits<SC>::magnitude(closestToRghtBoundDist)) {
667 rowSumDeviation = 0;
668 for (LocalOrdinal i = closestToLeftBound; i < closestToRghtBound ; i++) fixedSorted[i] = origSorted[i] - delta;
669 }
670 else {
671 rowSumDeviation = rowSumDeviation + closestToRghtBoundDist;
672// if (closestToRghtBound < nEntries) {
673 fixedSorted[closestToRghtBound] = origSorted[closestToRghtBound];
674 closestToRghtBound++;
675 // }
676 if (closestToRghtBound >= nEntries) closestToRghtBoundDist = aBigNumber;
677 else closestToRghtBoundDist = origSorted[closestToRghtBound] - rghtBound;
678 }
679 }
680 }
681
682 if (flipped) {
683 /* flip sign of fixedSorted and reverse ordering so that the positve version is sorted */
684
685 if ((nEntries%2) == 1) fixedSorted[(nEntries/2) ] = -fixedSorted[(nEntries/2) ];
686 for (LocalOrdinal i=0; i < nEntries/2; i++) {
687 temp=fixedSorted[i];
688 fixedSorted[i] = -fixedSorted[nEntries-1-i];
689 fixedSorted[nEntries-i-1] = -temp;
690 }
691 }
692 for (LocalOrdinal i = 0; i < nEntries; i++) fixedUnsorted[inds[i]] = fixedSorted[i];
693
694 /* check that no constraints are violated */
695
696 bool lowerViolation = false;
697 bool upperViolation = false;
698 bool sumViolation = false;
699 using TST = Teuchos::ScalarTraits<SC>;
700 temp = TST::zero();
701 for (LocalOrdinal i = 0; i < nEntries; i++) {
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];
705 }
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;
708
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??? ");
712
713 return hasFeasibleSol;
714}
715
716
717} //namespace MueLu
718
719#endif // MUELU_SAPFACTORY_DEF_HPP
720
721//TODO: restrictionMode_ should use the parameter list.
#define SET_VALID_ENTRY(name)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
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
Enforce constraints on prolongator.
void Build(Level &fineLevel, Level &coarseLevel) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void optimalSatisfyPConstraintsForScalarPDEs(RCP< Matrix > &P) const
bool constrainRow(Scalar *orig, LocalOrdinal nEntries, Scalar leftBound, Scalar rghtBound, Scalar rsumTarget, Scalar *fixedUnsorted, Scalar *scalarData) const
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetLumpedMatrixDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
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 > &params=Teuchos::null)
static Scalar PowerMethod(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixDiagonalInverse(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.
@ Runtime0
One-liner description of what is happening.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
scalar_type dampingFactor