MueLu Version of the Day
Loading...
Searching...
No Matches
BelosXpetraStatusTestGenResSubNorm.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MueLu: A package for multigrid based preconditioning
4//
5// Copyright 2012 NTESS and the MueLu contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP
11#define BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP
12
13#include "Xpetra_ConfigDefs.hpp"
14
15#include "Xpetra_BlockedCrsMatrix.hpp"
16
17#include "MueLu_Exceptions.hpp"
18
19#include <BelosConfigDefs.hpp>
20#include <BelosTypes.hpp>
21#include <BelosOperatorT.hpp>
22#include <BelosXpetraAdapterOperator.hpp>
23#include <BelosStatusTestGenResSubNorm.hpp>
24
25namespace Belos {
26
30template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
31class StatusTestGenResSubNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >, Teuchos::SerialDenseMatrix<LocalOrdinal, Scalar> >
32 : public StatusTestResNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >, Teuchos::SerialDenseMatrix<LocalOrdinal, Scalar> > {
33 public:
34 // Convenience typedefs
35 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> MV;
36 typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> BCRS;
37 typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> ME;
38 typedef Teuchos::SerialDenseMatrix<LocalOrdinal, Scalar> DM;
39 typedef Belos::OperatorT<MV> OP;
40
41 typedef Teuchos::ScalarTraits<Scalar> SCT;
42 typedef typename SCT::magnitudeType MagnitudeType;
43 typedef MultiVecTraits<Scalar, MV, DM> MVT;
44 typedef OperatorTraits<Scalar, MV, OP> OT;
45
47
48
61 StatusTestGenResSubNorm(MagnitudeType Tolerance, size_t subIdx, int quorum = -1, bool showMaxResNormOnly = false)
62 : tolerance_(Tolerance)
63 , subIdx_(subIdx)
64 , quorum_(quorum)
65 , showMaxResNormOnly_(showMaxResNormOnly)
66 , resnormtype_(TwoNorm)
67 , scaletype_(NormOfInitRes)
68 , scalenormtype_(TwoNorm)
69 , scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one())
70 , status_(Undefined)
71 , curBlksz_(0)
72 , curNumRHS_(0)
73 , curLSNum_(0)
74 , numrhs_(0)
75 , firstcallCheckStatus_(true)
76 , firstcallDefineResForm_(true)
77 , firstcallDefineScaleForm_(true)
78 , mapExtractor_(Teuchos::null) {}
79
83
85
86
88
94 int defineResForm(NormType TypeOfNorm) {
95 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_ == false, StatusTestError,
96 "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
97 firstcallDefineResForm_ = false;
98
99 resnormtype_ = TypeOfNorm;
100
101 return (0);
102 }
103
105
125 int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one()) {
126 TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_ == false, StatusTestError,
127 "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
128 firstcallDefineScaleForm_ = false;
129
130 scaletype_ = TypeOfScaling;
131 scalenormtype_ = TypeOfNorm;
132 scalevalue_ = ScaleValue;
133
134 return (0);
135 }
136
138
142 tolerance_ = tolerance;
143 return (0);
144 }
145
147
149 int setSubIdx(size_t subIdx) {
150 subIdx_ = subIdx;
151 return (0);
152 }
153
156 int setQuorum(int quorum) {
157 quorum_ = quorum;
158 return (0);
159 }
160
162 int setShowMaxResNormOnly(bool showMaxResNormOnly) {
163 showMaxResNormOnly_ = showMaxResNormOnly;
164 return (0);
165 }
166
168
170
171
178 StatusType checkStatus(Iteration<Scalar, MV, OP, DM>* iSolver) {
179 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
180 const LinearProblem<Scalar, MV, OP>& lp = iSolver->getProblem();
181 // Compute scaling term (done once for each block that's being solved)
182 if (firstcallCheckStatus_) {
183 StatusType status = firstCallCheckStatusSetup(iSolver);
184 if (status == Failed) {
185 status_ = Failed;
186 return (status_);
187 }
188 }
189
190 //
191 // This section computes the norm of the residual std::vector
192 //
193 if (curLSNum_ != lp.getLSNumber()) {
194 //
195 // We have moved on to the next rhs block
196 //
197 curLSNum_ = lp.getLSNumber();
198 curLSIdx_ = lp.getLSIndex();
199 curBlksz_ = (int)curLSIdx_.size();
200 int validLS = 0;
201 for (int i = 0; i < curBlksz_; ++i) {
202 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
203 validLS++;
204 }
205 curNumRHS_ = validLS;
206 curSoln_ = Teuchos::null;
207 //
208 } else {
209 //
210 // We are in the same rhs block, return if we are converged
211 //
212 if (status_ == Passed) {
213 return status_;
214 }
215 }
216
217 //
218 // Request the true residual for this block of right-hand sides.
219 //
220 Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
221 curSoln_ = lp.updateSolution(cur_update);
222 Teuchos::RCP<MV> cur_res = MVT::Clone(*curSoln_, MVT::GetNumberVecs(*curSoln_));
223 lp.computeCurrResVec(&*cur_res, &*curSoln_);
224 std::vector<MagnitudeType> tmp_resvector(MVT::GetNumberVecs(*cur_res));
225 MvSubNorm(*cur_res, subIdx_, tmp_resvector, resnormtype_);
226
227 typename std::vector<int>::iterator p = curLSIdx_.begin();
228 for (int i = 0; p < curLSIdx_.end(); ++p, ++i) {
229 // Check if this index is valid
230 if (*p != -1)
231 resvector_[*p] = tmp_resvector[i];
232 }
233
234 //
235 // Compute the new linear system residuals for testing.
236 // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
237 //
238 if (scalevector_.size() > 0) {
239 typename std::vector<int>::iterator pp = curLSIdx_.begin();
240 for (; pp < curLSIdx_.end(); ++pp) {
241 // Check if this index is valid
242 if (*pp != -1) {
243 // Scale the std::vector accordingly
244 if (scalevector_[*pp] != zero) {
245 // Don't intentionally divide by zero.
246 testvector_[*pp] = resvector_[*pp] / scalevector_[*pp] / scalevalue_;
247 } else {
248 testvector_[*pp] = resvector_[*pp] / scalevalue_;
249 }
250 }
251 }
252 } else {
253 typename std::vector<int>::iterator pp = curLSIdx_.begin();
254 for (; pp < curLSIdx_.end(); ++pp) {
255 // Check if this index is valid
256 if (*pp != -1)
257 testvector_[*pp] = resvector_[*pp] / scalevalue_;
258 }
259 }
260 // Check status of new linear system residuals and see if we have the quorum.
261 int have = 0;
262 ind_.resize(curLSIdx_.size());
263 typename std::vector<int>::iterator p2 = curLSIdx_.begin();
264 for (; p2 < curLSIdx_.end(); ++p2) {
265 // Check if this index is valid
266 if (*p2 != -1) {
267 // Check if any of the residuals are larger than the tolerance.
268 if (testvector_[*p2] > tolerance_) {
269 // do nothing.
270 } else if (testvector_[*p2] == Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero())) {
271 reset();
272 } else if (testvector_[*p2] <= tolerance_) {
273 ind_[have] = *p2;
274 have++;
275 } else {
276 // Throw an std::exception if a NaN is found.
277 status_ = Failed;
278 TEUCHOS_TEST_FOR_EXCEPTION(true, StatusTestError, "StatusTestGenResSubNorm::checkStatus(): NaN has been detected.");
279 }
280 }
281 }
282 ind_.resize(have);
283 int need = (quorum_ == -1) ? curNumRHS_ : quorum_;
284 status_ = (have >= need) ? Passed : Failed;
285 // Return the current status
286 return status_;
287 }
288
290 StatusType getStatus() const { return (status_); };
292
294
295
297 void reset() {
298 status_ = Undefined;
299 curBlksz_ = 0;
300 curLSNum_ = 0;
301 curLSIdx_.resize(0);
302 numrhs_ = 0;
303 ind_.resize(0);
304 firstcallCheckStatus_ = true;
305 curSoln_ = Teuchos::null;
306 }
307
309
311
312
314 void print(std::ostream& os, int indent = 0) const {
315 os.setf(std::ios_base::scientific);
316 for (int j = 0; j < indent; j++)
317 os << ' ';
318 printStatus(os, status_);
319 os << resFormStr();
320 if (status_ == Undefined)
321 os << ", tol = " << tolerance_ << std::endl;
322 else {
323 os << std::endl;
324 if (showMaxResNormOnly_ && curBlksz_ > 1) {
325 const MagnitudeType maxRelRes = *std::max_element(
326 testvector_.begin() + curLSIdx_[0], testvector_.begin() + curLSIdx_[curBlksz_ - 1]);
327 for (int j = 0; j < indent + 13; j++)
328 os << ' ';
329 os << "max{residual[" << curLSIdx_[0] << "..." << curLSIdx_[curBlksz_ - 1] << "]} = " << maxRelRes
330 << (maxRelRes <= tolerance_ ? " <= " : " > ") << tolerance_ << std::endl;
331 } else {
332 for (int i = 0; i < numrhs_; i++) {
333 for (int j = 0; j < indent + 13; j++)
334 os << ' ';
335 os << "residual [ " << i << " ] = " << testvector_[i];
336 os << ((testvector_[i] < tolerance_) ? " < " : (testvector_[i] == tolerance_) ? " == "
337 : (testvector_[i] > tolerance_) ? " > "
338 : " ")
339 << tolerance_ << std::endl;
340 }
341 }
342 }
343 os << std::endl;
344 }
345
347 void printStatus(std::ostream& os, StatusType type) const {
348 os << std::left << std::setw(13) << std::setfill('.');
349 switch (type) {
350 case Passed:
351 os << "Converged";
352 break;
353 case Failed:
354 os << "Unconverged";
355 break;
356 case Undefined:
357 default:
358 os << "**";
359 break;
360 }
361 os << std::left << std::setfill(' ');
362 return;
363 }
365
367
368
370 Teuchos::RCP<MV> getSolution() { return curSoln_; }
371
374 int getQuorum() const { return quorum_; }
375
377 size_t getSubIdx() const { return subIdx_; }
378
380 bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
381
383 std::vector<int> convIndices() { return ind_; }
384
386 MagnitudeType getTolerance() const { return (tolerance_); };
387
389 const std::vector<MagnitudeType>* getTestValue() const { return (&testvector_); };
390
392 const std::vector<MagnitudeType>* getResNormValue() const { return (&resvector_); };
393
395 const std::vector<MagnitudeType>* getScaledNormValue() const { return (&scalevector_); };
396
399 bool getLOADetected() const { return false; }
400
402
405
411 StatusType firstCallCheckStatusSetup(Iteration<Scalar, MV, OP>* iSolver) {
412 int i;
413 MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
414 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
415 const LinearProblem<Scalar, MV, OP>& lp = iSolver->getProblem();
416 // Compute scaling term (done once for each block that's being solved)
417 if (firstcallCheckStatus_) {
418 //
419 // Get some current solver information.
420 //
421 firstcallCheckStatus_ = false;
422
423 // try to access the underlying blocked operator
424 Teuchos::RCP<const OP> Op = lp.getOperator();
425 Teuchos::RCP<const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xOp =
426 Teuchos::rcp_dynamic_cast<const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(Op);
427 TEUCHOS_TEST_FOR_EXCEPTION(xOp.is_null(), MueLu::Exceptions::BadCast, "Bad cast from \'const Belos::OperatorT\' to \'const Belos::XpetraOp\'. The origin type is " << typeid(const OP).name() << ".");
428 Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xIntOp =
429 xOp->getOperator();
430 TEUCHOS_TEST_FOR_EXCEPTION(xIntOp.is_null(), MueLu::Exceptions::BadCast, "Cannot access Xpetra::Operator stored in Belos::XpetraOperator.");
431 Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xMat =
432 Teuchos::rcp_dynamic_cast<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xIntOp);
433 TEUCHOS_TEST_FOR_EXCEPTION(xMat.is_null(), MueLu::Exceptions::RuntimeError, "Cannot access Xpetra::Matrix stored in Belos::XpetraOp. Error.");
434 Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bMat = Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xMat);
435 TEUCHOS_TEST_FOR_EXCEPTION(bMat.is_null(), MueLu::Exceptions::BadCast, "Bad cast from \'const Xpetra::Matrix\' to \'const Xpetra::BlockedCrsMatrix\'. The origin type is " << typeid(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>).name() << ". Note: you need a BlockedCrsMatrix object for the StatusTestGenResSubNorm to work!");
436 mapExtractor_ = bMat->getRangeMapExtractor();
437 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_.is_null(), MueLu::Exceptions::RuntimeError, "Could not extract map extractor from BlockedCrsMatrix. Error.");
438 TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_->NumMaps() <= subIdx_, MueLu::Exceptions::RuntimeError, "The multivector is only split into " << mapExtractor_->NumMaps() << " sub parts. Cannot access sub-block " << subIdx_ << ".");
439
440 // calculate initial norms
441 if (scaletype_ == NormOfRHS) {
442 Teuchos::RCP<const MV> rhs = lp.getRHS();
443 numrhs_ = MVT::GetNumberVecs(*rhs);
444 scalevector_.resize(numrhs_);
445 MvSubNorm(*rhs, subIdx_, scalevector_, scalenormtype_);
446 } else if (scaletype_ == NormOfInitRes) {
447 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
448 numrhs_ = MVT::GetNumberVecs(*init_res);
449 scalevector_.resize(numrhs_);
450 MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
451 } else if (scaletype_ == NormOfPrecInitRes) {
452 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
453 numrhs_ = MVT::GetNumberVecs(*init_res);
454 scalevector_.resize(numrhs_);
455 MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
456 } else if (scaletype_ == NormOfFullInitRes) {
457 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
458 numrhs_ = MVT::GetNumberVecs(*init_res);
459 scalevector_.resize(numrhs_);
460 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
461 scalevalue_ = one;
462 } else if (scaletype_ == NormOfFullPrecInitRes) {
463 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
464 numrhs_ = MVT::GetNumberVecs(*init_res);
465 scalevector_.resize(numrhs_);
466 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
467 scalevalue_ = one;
468 } else if (scaletype_ == NormOfFullScaledInitRes) {
469 Teuchos::RCP<const MV> init_res = lp.getInitResVec();
470 numrhs_ = MVT::GetNumberVecs(*init_res);
471 scalevector_.resize(numrhs_);
472 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
473 MvScalingRatio(*init_res, subIdx_, scalevalue_);
474 } else if (scaletype_ == NormOfFullScaledPrecInitRes) {
475 Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
476 numrhs_ = MVT::GetNumberVecs(*init_res);
477 scalevector_.resize(numrhs_);
478 MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
479 MvScalingRatio(*init_res, subIdx_, scalevalue_);
480 } else {
481 numrhs_ = MVT::GetNumberVecs(*(lp.getRHS()));
482 }
483
484 resvector_.resize(numrhs_);
485 testvector_.resize(numrhs_);
486
487 curLSNum_ = lp.getLSNumber();
488 curLSIdx_ = lp.getLSIndex();
489 curBlksz_ = (int)curLSIdx_.size();
490 int validLS = 0;
491 for (i = 0; i < curBlksz_; ++i) {
492 if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
493 validLS++;
494 }
495 curNumRHS_ = validLS;
496 //
497 // Initialize the testvector.
498 for (i = 0; i < numrhs_; i++) {
499 testvector_[i] = one;
500 }
501
502 // Return an error if the scaling is zero.
503 if (scalevalue_ == zero) {
504 return Failed;
505 }
506 }
507 return Undefined;
508 }
510
513
515 std::string description() const {
516 std::ostringstream oss;
517 oss << "Belos::StatusTestGenResSubNorm<>: " << resFormStr();
518 oss << ", tol = " << tolerance_;
519 return oss.str();
520 }
522
523 protected:
524 private:
526
527
528 std::string resFormStr() const {
529 std::ostringstream oss;
530 oss << "(";
531 oss << ((resnormtype_ == OneNorm) ? "1-Norm" : (resnormtype_ == TwoNorm) ? "2-Norm"
532 : "Inf-Norm");
533 oss << " Exp";
534 oss << " Res Vec [" << subIdx_ << "]) ";
535
536 // If there is no residual scaling, return current string.
537 if (scaletype_ != None) {
538 // Insert division sign.
539 oss << "/ ";
540
541 // Determine output string for scaling, if there is any.
542 if (scaletype_ == UserProvided)
543 oss << " (User Scale)";
544 else {
545 oss << "(";
546 oss << ((scalenormtype_ == OneNorm) ? "1-Norm" : (resnormtype_ == TwoNorm) ? "2-Norm"
547 : "Inf-Norm");
548 if (scaletype_ == NormOfInitRes)
549 oss << " Res0 [" << subIdx_ << "]";
550 else if (scaletype_ == NormOfPrecInitRes)
551 oss << " Prec Res0 [" << subIdx_ << "]";
552 else if (scaletype_ == NormOfFullInitRes)
553 oss << " Full Res0 [" << subIdx_ << "]";
554 else if (scaletype_ == NormOfFullPrecInitRes)
555 oss << " Full Prec Res0 [" << subIdx_ << "]";
556 else if (scaletype_ == NormOfFullScaledInitRes)
557 oss << " scaled Full Res0 [" << subIdx_ << "]";
558 else if (scaletype_ == NormOfFullScaledPrecInitRes)
559 oss << " scaled Full Prec Res0 [" << subIdx_ << "]";
560 else
561 oss << " RHS [" << subIdx_ << "]";
562 oss << ")";
563 }
564 }
565
566 // TODO add a tagging name
567
568 return oss.str();
569 }
570
572
574
575
576 // calculate norm of partial multivector
577 void MvSubNorm(const MV& mv, size_t block, std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& normVec, NormType type = TwoNorm) {
578 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
579
580 Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
581 MVT::MvNorm(*SubVec, normVec, type);
582 }
583
584 // calculate ration of sub-vector length to full vector length (for scalevalue_)
585 void MvScalingRatio(const MV& mv, size_t block, MagnitudeType& lengthRatio) {
586 Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
587
588 Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
589
590 lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
591 }
593
595
596
599
601 size_t subIdx_;
602
605
608
610 NormType resnormtype_;
611
613 ScaleType scaletype_;
614
617
620
622 std::vector<MagnitudeType> scalevector_;
623
625 std::vector<MagnitudeType> resvector_;
626
628 std::vector<MagnitudeType> testvector_;
629
631 std::vector<int> ind_;
632
634 Teuchos::RCP<MV> curSoln_;
635
637 StatusType status_;
638
641
644
646 std::vector<int> curLSIdx_;
647
650
653
656
659
662
664 Teuchos::RCP<const ME> mapExtractor_;
666};
667
668} // namespace Belos
669
670#endif /* BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP */
MueLu::DefaultScalar Scalar
int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue=Teuchos::ScalarTraits< MagnitudeType >::one())
Define form of the scaling, its norm, its optional weighting std::vector, or, alternatively,...
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
magnitude_type tolerance