waves
Basic FE playground
BelosXpetraStatusTestGenResSubNorm_MP_Vector.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 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_MP_VECTORHPP
48 #define BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_MP_VECTORHPP
49 
50 #include "Stokhos_config.h"
51 
52 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
53 
54 #include "Xpetra_ConfigDefs.hpp"
55 
56 #include "Xpetra_BlockedCrsMatrix.hpp"
57 
58 #include "MueLu_Exceptions.hpp"
59 
60 #include <BelosConfigDefs.hpp>
61 #include <BelosTypes.hpp>
62 #include <BelosOperatorT.hpp>
63 #include <BelosXpetraAdapterOperator.hpp>
64 #include <BelosStatusTestGenResSubNorm.hpp>
65 #include <BelosXpetraStatusTestGenResSubNorm.hpp>
66 #include "EnsembleTraits.h"
67 
68 namespace Belos
69 {
70 
74 template <class Storage, class LocalOrdinal, class GlobalOrdinal, class Node>
75 class StatusTestGenResSubNorm<typename Teuchos::ScalarTraits<Sacado::MP::Vector<Storage>>::magnitudeType, Xpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal, GlobalOrdinal, Node>>>
76  : public StatusTestResNorm<typename Teuchos::ScalarTraits<Sacado::MP::Vector<Storage>>::magnitudeType, Xpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal, GlobalOrdinal, Node>>>
77 {
78 
79 public:
80  // Convenience typedefs
81  typedef Sacado::MP::Vector<Storage> ScalarType;
82 
83  typedef Xpetra::MultiVector<ScalarType, LocalOrdinal, GlobalOrdinal, Node> MV;
84  typedef Xpetra::BlockedCrsMatrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node> BCRS;
85  typedef Xpetra::MapExtractor<ScalarType, LocalOrdinal, GlobalOrdinal, Node> ME;
86  typedef Belos::OperatorT<MV> OP;
87 
88  typedef Teuchos::ScalarTraits<ScalarType> SCT;
89  typedef typename SCT::magnitudeType MagnitudeType;
90  typedef MultiVecTraits<ScalarType, MV> MVT;
91  typedef OperatorTraits<ScalarType, MV, OP> OT;
92 
94 
95 
108  StatusTestGenResSubNorm(MagnitudeType Tolerance, size_t subIdx, int quorum = -1, bool showMaxResNormOnly = false)
109  : tolerance_(Tolerance),
110  subIdx_(subIdx),
111  quorum_(quorum),
112  showMaxResNormOnly_(showMaxResNormOnly),
113  resnormtype_(TwoNorm),
114  scaletype_(NormOfInitRes),
115  scalenormtype_(TwoNorm),
116  scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one()),
117  status_(Undefined),
118  curBlksz_(0),
119  curNumRHS_(0),
120  curLSNum_(0),
121  numrhs_(0),
122  firstcallCheckStatus_(true),
123  firstcallDefineResForm_(true),
124  firstcallDefineScaleForm_(true),
125  mapExtractor_(Teuchos::null),
126  ensemble_iterations(EnsembleTraits<MagnitudeType>::size, 0) {}
127 
129  virtual ~StatusTestGenResSubNorm(){};
131 
133 
134 
136 
142  int defineResForm(NormType TypeOfNorm)
143  {
144  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_ == false, StatusTestError,
145  "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
146  firstcallDefineResForm_ = false;
147 
148  resnormtype_ = TypeOfNorm;
149 
150  return (0);
151  }
152 
154 
174  int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one())
175  {
176  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_ == false, StatusTestError,
177  "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
178  firstcallDefineScaleForm_ = false;
179 
180  scaletype_ = TypeOfScaling;
181  scalenormtype_ = TypeOfNorm;
182  scalevalue_ = ScaleValue;
183 
184  return (0);
185  }
186 
188 
191  int setTolerance(MagnitudeType tolerance)
192  {
193  tolerance_ = tolerance;
194  return (0);
195  }
196 
198 
200  int setSubIdx(size_t subIdx)
201  {
202  subIdx_ = subIdx;
203  return (0);
204  }
205 
208  int setQuorum(int quorum)
209  {
210  quorum_ = quorum;
211  return (0);
212  }
213 
215  int setShowMaxResNormOnly(bool showMaxResNormOnly)
216  {
217  showMaxResNormOnly_ = showMaxResNormOnly;
218  return (0);
219  }
220 
222 
224 
225 
232  StatusType checkStatus(Iteration<MagnitudeType, MV, OP> *iSolver)
233  {
235  const int ensemble_size = ET::size;
236  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
237  const LinearProblem<MagnitudeType, MV, OP> &lp = iSolver->getProblem();
238  // Compute scaling term (done once for each block that's being solved)
239  if (firstcallCheckStatus_)
240  {
241  StatusType status = firstCallCheckStatusSetup(iSolver);
242  if (status == Failed)
243  {
244  status_ = Failed;
245  return (status_);
246  }
247  }
248 
249  //
250  // This section computes the norm of the residual std::vector
251  //
252  if (curLSNum_ != lp.getLSNumber())
253  {
254  //
255  // We have moved on to the next rhs block
256  //
257  curLSNum_ = lp.getLSNumber();
258  curLSIdx_ = lp.getLSIndex();
259  curBlksz_ = (int)curLSIdx_.size();
260  int validLS = 0;
261  for (int i = 0; i < curBlksz_; ++i)
262  {
263  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
264  validLS++;
265  }
266  curNumRHS_ = validLS;
267  curSoln_ = Teuchos::null;
268  //
269  }
270  else
271  {
272  //
273  // We are in the same rhs block, return if we are converged
274  //
275  if (status_ == Passed)
276  {
277  return status_;
278  }
279  }
280 
281  //
282  // Request the true residual for this block of right-hand sides.
283  //
284  Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
285  curSoln_ = lp.updateSolution(cur_update);
286  Teuchos::RCP<MV> cur_res = MVT::Clone(*curSoln_, MVT::GetNumberVecs(*curSoln_));
287  lp.computeCurrResVec(&*cur_res, &*curSoln_);
288  std::vector<MagnitudeType> tmp_resvector(MVT::GetNumberVecs(*cur_res));
289  MvSubNorm(*cur_res, subIdx_, tmp_resvector, resnormtype_);
290 
291  typename std::vector<int>::iterator p = curLSIdx_.begin();
292  for (int i = 0; p < curLSIdx_.end(); ++p, ++i)
293  {
294  // Check if this index is valid
295  if (*p != -1)
296  resvector_[*p] = tmp_resvector[i];
297  }
298 
299  //
300  // Compute the new linear system residuals for testing.
301  // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
302  //
303  if (scalevector_.size() > 0)
304  {
305  typename std::vector<int>::iterator pp = curLSIdx_.begin();
306  for (; pp < curLSIdx_.end(); ++pp)
307  {
308  // Check if this index is valid
309  if (*pp != -1)
310  {
311  // Scale the std::vector accordingly
312  if (scalevector_[*pp] != zero)
313  {
314  // Don't intentionally divide by zero.
315  testvector_[*pp] = resvector_[*pp] / scalevector_[*pp] / scalevalue_;
316  }
317  else
318  {
319  testvector_[*pp] = resvector_[*pp] / scalevalue_;
320  }
321  }
322  }
323  }
324  else
325  {
326  typename std::vector<int>::iterator pp = curLSIdx_.begin();
327  for (; pp < curLSIdx_.end(); ++pp)
328  {
329  // Check if this index is valid
330  if (*pp != -1)
331  testvector_[*pp] = resvector_[*pp] / scalevalue_;
332  }
333  }
334  // Check status of new linear system residuals and see if we have the quorum.
335  int have = 0;
336  ind_.resize(curLSIdx_.size());
337  typename std::vector<int>::iterator p2 = curLSIdx_.begin();
338  for (; p2 < curLSIdx_.end(); ++p2)
339  {
340  // Check if this index is valid
341  if (*p2 != -1)
342  {
343  // Check if any of the residuals are larger than the tolerance.
344  bool all_converged = true;
345  for (int ii = 0; ii < ensemble_size; ++ii)
346  {
347  if (ET::coeff(testvector_[*p2], ii) > ET::coeff(tolerance_, ii))
348  {
349  ++ensemble_iterations[ii];
350  all_converged = false;
351  }
352  else if (!(ET::coeff(testvector_[*p2], ii) <= ET::coeff(tolerance_, ii)))
353  {
354  // Throw an std::exception if the current residual norm is
355  // NaN. We know that it's NaN because it is not less than,
356  // equal to, or greater than the current tolerance. This is
357  // only possible if either the residual norm or the current
358  // tolerance is NaN; we assume the former. We also mark the
359  // test as failed, in case you want to catch the exception.
360  status_ = Failed;
361  TEUCHOS_TEST_FOR_EXCEPTION(true, StatusTestError, "StatusTestGenResSubNorm::checkStatus(): NaN has been detected.");
362  }
363  }
364  if (all_converged)
365  {
366  ind_[have] = *p2;
367  have++;
368  }
369  }
370  }
371  ind_.resize(have);
372  int need = (quorum_ == -1) ? curNumRHS_ : quorum_;
373  status_ = (have >= need) ? Passed : Failed;
374  // Return the current status
375  return status_;
376  }
377 
379  StatusType getStatus() const { return (status_); };
381 
383 
384 
386  void reset()
387  {
388  status_ = Undefined;
389  curBlksz_ = 0;
390  curLSNum_ = 0;
391  curLSIdx_.resize(0);
392  numrhs_ = 0;
393  ind_.resize(0);
394  firstcallCheckStatus_ = true;
395  curSoln_ = Teuchos::null;
396  }
397 
399 
401 
402 
404  void print(std::ostream &os, int indent = 0) const
405  {
406  os.setf(std::ios_base::scientific);
407  for (int j = 0; j < indent; j++)
408  os << ' ';
409  printStatus(os, status_);
410  os << resFormStr();
411  if (status_ == Undefined)
412  os << ", tol = " << tolerance_ << std::endl;
413  else
414  {
415  os << std::endl;
416  if (showMaxResNormOnly_ && curBlksz_ > 1)
417  {
418  const MagnitudeType maxRelRes = *std::max_element(
419  testvector_.begin() + curLSIdx_[0], testvector_.begin() + curLSIdx_[curBlksz_ - 1]);
420  for (int j = 0; j < indent + 13; j++)
421  os << ' ';
422  os << "max{residual[" << curLSIdx_[0] << "..." << curLSIdx_[curBlksz_ - 1] << "]} = " << maxRelRes
423  << (maxRelRes <= tolerance_ ? " <= " : " > ") << tolerance_ << std::endl;
424  }
425  else
426  {
427  for (int i = 0; i < numrhs_; i++)
428  {
429  for (int j = 0; j < indent + 13; j++)
430  os << ' ';
431  os << "residual [ " << i << " ] = " << testvector_[i];
432  os << ((testvector_[i] < tolerance_) ? " < " : (testvector_[i] == tolerance_) ? " == " : (testvector_[i] > tolerance_) ? " > " : " ") << tolerance_ << std::endl;
433  }
434  }
435  }
436  os << std::endl;
437  }
438 
440  void printStatus(std::ostream &os, StatusType type) const
441  {
442  os << std::left << std::setw(13) << std::setfill('.');
443  switch (type)
444  {
445  case Passed:
446  os << "Converged";
447  break;
448  case Failed:
449  os << "Unconverged";
450  break;
451  case Undefined:
452  default:
453  os << "**";
454  break;
455  }
456  os << std::left << std::setfill(' ');
457  return;
458  }
460 
462 
463 
465  Teuchos::RCP<MV> getSolution() { return curSoln_; }
466 
469  int getQuorum() const { return quorum_; }
470 
472  size_t getSubIdx() const { return subIdx_; }
473 
475  bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
476 
478  std::vector<int> convIndices() { return ind_; }
479 
481  MagnitudeType getTolerance() const { return (tolerance_); };
482 
484  const std::vector<MagnitudeType> *getTestValue() const { return (&testvector_); };
485 
487  const std::vector<MagnitudeType> *getResNormValue() const { return (&resvector_); };
488 
490  const std::vector<MagnitudeType> *getScaledNormValue() const { return (&scalevector_); };
491 
494  bool getLOADetected() const { return false; }
495 
496  const std::vector<int> getEnsembleIterations() const { return ensemble_iterations; }
497 
499 
502 
508  StatusType firstCallCheckStatusSetup(Iteration<MagnitudeType, MV, OP> *iSolver)
509  {
510  int i;
511  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
512  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
513  const LinearProblem<MagnitudeType, MV, OP> &lp = iSolver->getProblem();
514  // Compute scaling term (done once for each block that's being solved)
515  if (firstcallCheckStatus_)
516  {
517  //
518  // Get some current solver information.
519  //
520  firstcallCheckStatus_ = false;
521 
522  // try to access the underlying blocked operator
523  Teuchos::RCP<const OP> Op = lp.getOperator();
524  Teuchos::RCP<const Belos::XpetraOp<ScalarType, LocalOrdinal, GlobalOrdinal, Node>> xOp =
525  Teuchos::rcp_dynamic_cast<const Belos::XpetraOp<ScalarType, LocalOrdinal, GlobalOrdinal, Node>>(Op);
526  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() << ".");
527  Teuchos::RCP<const Xpetra::Operator<ScalarType, LocalOrdinal, GlobalOrdinal, Node>> xIntOp =
528  xOp->getOperator();
529  TEUCHOS_TEST_FOR_EXCEPTION(xIntOp.is_null(), MueLu::Exceptions::BadCast, "Cannot access Xpetra::Operator stored in Belos::XpetraOperator.");
530  Teuchos::RCP<const Xpetra::Matrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>> xMat =
531  Teuchos::rcp_dynamic_cast<const Xpetra::Matrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>>(xIntOp);
532  TEUCHOS_TEST_FOR_EXCEPTION(xMat.is_null(), MueLu::Exceptions::RuntimeError, "Cannot access Xpetra::Matrix stored in Belos::XpetraOp. Error.");
533  Teuchos::RCP<const Xpetra::BlockedCrsMatrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>> bMat = Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<ScalarType, LocalOrdinal, GlobalOrdinal, Node>>(xMat);
534  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<ScalarType, LocalOrdinal, GlobalOrdinal, Node>).name() << ". Note: you need a BlockedCrsMatrix object for the StatusTestGenResSubNorm to work!");
535  mapExtractor_ = bMat->getRangeMapExtractor();
536  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_.is_null(), MueLu::Exceptions::RuntimeError, "Could not extract map extractor from BlockedCrsMatrix. Error.");
537  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_ << ".");
538 
539  // calculate initial norms
540  if (scaletype_ == NormOfRHS)
541  {
542  Teuchos::RCP<const MV> rhs = lp.getRHS();
543  numrhs_ = MVT::GetNumberVecs(*rhs);
544  scalevector_.resize(numrhs_);
545  MvSubNorm(*rhs, subIdx_, scalevector_, scalenormtype_);
546  }
547  else if (scaletype_ == NormOfInitRes)
548  {
549  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
550  numrhs_ = MVT::GetNumberVecs(*init_res);
551  scalevector_.resize(numrhs_);
552  MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
553  }
554  else if (scaletype_ == NormOfPrecInitRes)
555  {
556  Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
557  numrhs_ = MVT::GetNumberVecs(*init_res);
558  scalevector_.resize(numrhs_);
559  MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
560  }
561  else if (scaletype_ == NormOfFullInitRes)
562  {
563  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
564  numrhs_ = MVT::GetNumberVecs(*init_res);
565  scalevector_.resize(numrhs_);
566  MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
567  scalevalue_ = one;
568  }
569  else if (scaletype_ == NormOfFullPrecInitRes)
570  {
571  Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
572  numrhs_ = MVT::GetNumberVecs(*init_res);
573  scalevector_.resize(numrhs_);
574  MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
575  scalevalue_ = one;
576  }
577  else if (scaletype_ == NormOfFullScaledInitRes)
578  {
579  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
580  numrhs_ = MVT::GetNumberVecs(*init_res);
581  scalevector_.resize(numrhs_);
582  MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
583  MvScalingRatio(*init_res, subIdx_, scalevalue_);
584  }
585  else if (scaletype_ == NormOfFullScaledPrecInitRes)
586  {
587  Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
588  numrhs_ = MVT::GetNumberVecs(*init_res);
589  scalevector_.resize(numrhs_);
590  MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
591  MvScalingRatio(*init_res, subIdx_, scalevalue_);
592  }
593  else
594  {
595  numrhs_ = MVT::GetNumberVecs(*(lp.getRHS()));
596  }
597 
598  resvector_.resize(numrhs_);
599  testvector_.resize(numrhs_);
600 
601  curLSNum_ = lp.getLSNumber();
602  curLSIdx_ = lp.getLSIndex();
603  curBlksz_ = (int)curLSIdx_.size();
604  int validLS = 0;
605  for (i = 0; i < curBlksz_; ++i)
606  {
607  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
608  validLS++;
609  }
610  curNumRHS_ = validLS;
611  //
612  // Initialize the testvector.
613  for (i = 0; i < numrhs_; i++)
614  {
615  testvector_[i] = one;
616  }
617 
618  // Return an error if the scaling is zero.
619  if (scalevalue_ == zero)
620  {
621  return Failed;
622  }
623  }
624  return Undefined;
625  }
627 
630 
632  std::string description() const
633  {
634  std::ostringstream oss;
635  oss << "Belos::StatusTestGenResSubNorm<>: " << resFormStr();
636  oss << ", tol = " << tolerance_;
637  return oss.str();
638  }
640 
641 protected:
642 private:
644 
645 
646  std::string resFormStr() const
647  {
648  std::ostringstream oss;
649  oss << "(";
650  oss << ((resnormtype_ == OneNorm) ? "1-Norm" : (resnormtype_ == TwoNorm) ? "2-Norm" : "Inf-Norm");
651  oss << " Exp";
652  oss << " Res Vec [" << subIdx_ << "]) ";
653 
654  // If there is no residual scaling, return current string.
655  if (scaletype_ != None)
656  {
657  // Insert division sign.
658  oss << "/ ";
659 
660  // Determine output string for scaling, if there is any.
661  if (scaletype_ == UserProvided)
662  oss << " (User Scale)";
663  else
664  {
665  oss << "(";
666  oss << ((scalenormtype_ == OneNorm) ? "1-Norm" : (resnormtype_ == TwoNorm) ? "2-Norm" : "Inf-Norm");
667  if (scaletype_ == NormOfInitRes)
668  oss << " Res0 [" << subIdx_ << "]";
669  else if (scaletype_ == NormOfPrecInitRes)
670  oss << " Prec Res0 [" << subIdx_ << "]";
671  else if (scaletype_ == NormOfFullInitRes)
672  oss << " Full Res0 [" << subIdx_ << "]";
673  else if (scaletype_ == NormOfFullPrecInitRes)
674  oss << " Full Prec Res0 [" << subIdx_ << "]";
675  else if (scaletype_ == NormOfFullScaledInitRes)
676  oss << " scaled Full Res0 [" << subIdx_ << "]";
677  else if (scaletype_ == NormOfFullScaledPrecInitRes)
678  oss << " scaled Full Prec Res0 [" << subIdx_ << "]";
679  else
680  oss << " RHS [" << subIdx_ << "]";
681  oss << ")";
682  }
683  }
684 
685  // TODO add a tagging name
686 
687  return oss.str();
688  }
689 
691 
693 
694 
695  // calculate norm of partial multivector
696  void MvSubNorm(const MV &mv, size_t block, std::vector<MagnitudeType> &normVec, NormType type = TwoNorm)
697  {
698 
699  Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
700 
701  Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
702  typedef MultiVecTraits<ScalarType, MV> MVT;
703  MVT::MvNorm(*SubVec, normVec, type);
704  }
705 
706  // calculate ration of sub-vector length to full vector length (for scalevalue_)
707  void MvScalingRatio(const MV &mv, size_t block, MagnitudeType &lengthRatio)
708  {
709  Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
710 
711  Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
712 
713  lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
714  }
716 
718 
719 
721  MagnitudeType tolerance_;
722 
724  size_t subIdx_;
725 
727  int quorum_;
728 
730  bool showMaxResNormOnly_;
731 
733  NormType resnormtype_;
734 
736  ScaleType scaletype_;
737 
739  NormType scalenormtype_;
740 
742  MagnitudeType scalevalue_;
743 
745  std::vector<MagnitudeType> scalevector_;
746 
748  std::vector<MagnitudeType> resvector_;
749 
751  std::vector<MagnitudeType> testvector_;
752 
754  std::vector<int> ind_;
755 
757  Teuchos::RCP<MV> curSoln_;
758 
760  StatusType status_;
761 
763  int curBlksz_;
764 
766  int curNumRHS_;
767 
769  std::vector<int> curLSIdx_;
770 
772  int curLSNum_;
773 
775  int numrhs_;
776 
778  bool firstcallCheckStatus_;
779 
781  bool firstcallDefineResForm_;
782 
784  bool firstcallDefineScaleForm_;
785 
787  Teuchos::RCP<const ME> mapExtractor_;
788 
790  std::vector<int> ensemble_iterations;
791 
793 };
794 
795 } // namespace Belos
796 
797 #else
798 
799 #include "Xpetra_ConfigDefs.hpp"
800 
801 #include "Xpetra_BlockedCrsMatrix.hpp"
802 
803 #include "MueLu_Exceptions.hpp"
804 
805 #include <BelosConfigDefs.hpp>
806 #include <BelosTypes.hpp>
807 #include <BelosOperatorT.hpp>
808 #include <BelosXpetraAdapterOperator.hpp>
809 #include <BelosStatusTestGenResSubNorm.hpp>
810 #include "EnsembleTraits.h"
811 
812 namespace Belos
813 {
814 
818 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
819 class StatusTestGenResSubNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>>
820  : public StatusTestResNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>>
821 {
822 
823 public:
824  // Convenience typedefs
825  typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> MV;
826  typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> BCRS;
827  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> ME;
828  typedef Belos::OperatorT<MV> OP;
829 
830  typedef Teuchos::ScalarTraits<Scalar> SCT;
831  typedef typename SCT::magnitudeType MagnitudeType;
832  typedef MultiVecTraits<Scalar, MV> MVT;
833  typedef OperatorTraits<Scalar, MV, OP> OT;
834 
836 
837 
848  StatusTestGenResSubNorm(MagnitudeType Tolerance, size_t subIdx, int quorum = -1, bool showMaxResNormOnly = false)
849  : tolerance_(Tolerance),
850  subIdx_(subIdx),
851  quorum_(quorum),
852  showMaxResNormOnly_(showMaxResNormOnly),
853  resnormtype_(TwoNorm),
854  scaletype_(NormOfInitRes),
855  scalenormtype_(TwoNorm),
856  scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one()),
857  status_(Undefined),
858  curBlksz_(0),
859  curNumRHS_(0),
860  curLSNum_(0),
861  numrhs_(0),
862  firstcallCheckStatus_(true),
863  firstcallDefineResForm_(true),
864  firstcallDefineScaleForm_(true),
865  mapExtractor_(Teuchos::null),
866  ensemble_iterations(EnsembleTraits<Scalar>::size, 0) {}
867 
869  virtual ~StatusTestGenResSubNorm(){};
871 
873 
874 
876 
882  int defineResForm(NormType TypeOfNorm)
883  {
884  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_ == false, StatusTestError,
885  "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
886  firstcallDefineResForm_ = false;
887 
888  resnormtype_ = TypeOfNorm;
889 
890  return (0);
891  }
892 
894 
913  int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one())
914  {
915  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_ == false, StatusTestError,
916  "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
917  firstcallDefineScaleForm_ = false;
918 
919  scaletype_ = TypeOfScaling;
920  scalenormtype_ = TypeOfNorm;
921  scalevalue_ = ScaleValue;
922 
923  return (0);
924  }
925 
927 
930  int setTolerance(MagnitudeType tolerance)
931  {
932  tolerance_ = tolerance;
933  return (0);
934  }
935 
937 
939  int setSubIdx(size_t subIdx)
940  {
941  subIdx_ = subIdx;
942  return (0);
943  }
944 
947  int setQuorum(int quorum)
948  {
949  quorum_ = quorum;
950  return (0);
951  }
952 
954  int setShowMaxResNormOnly(bool showMaxResNormOnly)
955  {
956  showMaxResNormOnly_ = showMaxResNormOnly;
957  return (0);
958  }
959 
961 
963 
964 
970  StatusType checkStatus(Iteration<Scalar, MV, OP> *iSolver)
971  {
972  typedef EnsembleTraits<Scalar> ET;
973  const int ensemble_size = ET::size;
974  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
975  const LinearProblem<Scalar, MV, OP> &lp = iSolver->getProblem();
976  // Compute scaling term (done once for each block that's being solved)
977  if (firstcallCheckStatus_)
978  {
979  StatusType status = firstCallCheckStatusSetup(iSolver);
980  if (status == Failed)
981  {
982  status_ = Failed;
983  return (status_);
984  }
985  }
986 
987  //
988  // This section computes the norm of the residual std::vector
989  //
990  if (curLSNum_ != lp.getLSNumber())
991  {
992  //
993  // We have moved on to the next rhs block
994  //
995  curLSNum_ = lp.getLSNumber();
996  curLSIdx_ = lp.getLSIndex();
997  curBlksz_ = (int)curLSIdx_.size();
998  int validLS = 0;
999  for (int i = 0; i < curBlksz_; ++i)
1000  {
1001  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
1002  validLS++;
1003  }
1004  curNumRHS_ = validLS;
1005  curSoln_ = Teuchos::null;
1006  //
1007  }
1008  else
1009  {
1010  //
1011  // We are in the same rhs block, return if we are converged
1012  //
1013  if (status_ == Passed)
1014  {
1015  return status_;
1016  }
1017  }
1018 
1019  //
1020  // Request the true residual for this block of right-hand sides.
1021  //
1022  Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
1023  curSoln_ = lp.updateSolution(cur_update);
1024  Teuchos::RCP<MV> cur_res = MVT::Clone(*curSoln_, MVT::GetNumberVecs(*curSoln_));
1025  lp.computeCurrResVec(&*cur_res, &*curSoln_);
1026  std::vector<MagnitudeType> tmp_resvector(MVT::GetNumberVecs(*cur_res));
1027  MvSubNorm(*cur_res, subIdx_, tmp_resvector, resnormtype_);
1028 
1029  typename std::vector<int>::iterator p = curLSIdx_.begin();
1030  for (int i = 0; p < curLSIdx_.end(); ++p, ++i)
1031  {
1032  // Check if this index is valid
1033  if (*p != -1)
1034  resvector_[*p] = tmp_resvector[i];
1035  }
1036 
1037  //
1038  // Compute the new linear system residuals for testing.
1039  // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
1040  //
1041  if (scalevector_.size() > 0)
1042  {
1043  typename std::vector<int>::iterator pp = curLSIdx_.begin();
1044  for (; pp < curLSIdx_.end(); ++pp)
1045  {
1046  // Check if this index is valid
1047  if (*pp != -1)
1048  {
1049  // Scale the std::vector accordingly
1050  if (scalevector_[*pp] != zero)
1051  {
1052  // Don't intentionally divide by zero.
1053  testvector_[*pp] = resvector_[*pp] / scalevector_[*pp] / scalevalue_;
1054  }
1055  else
1056  {
1057  testvector_[*pp] = resvector_[*pp] / scalevalue_;
1058  }
1059  }
1060  }
1061  }
1062  else
1063  {
1064  typename std::vector<int>::iterator pp = curLSIdx_.begin();
1065  for (; pp < curLSIdx_.end(); ++pp)
1066  {
1067  // Check if this index is valid
1068  if (*pp != -1)
1069  testvector_[*pp] = resvector_[*pp] / scalevalue_;
1070  }
1071  }
1072  // Check status of new linear system residuals and see if we have the quorum.
1073  int have = 0;
1074  ind_.resize(curLSIdx_.size());
1075  typename std::vector<int>::iterator p2 = curLSIdx_.begin();
1076  for (; p2 < curLSIdx_.end(); ++p2)
1077  {
1078  // Check if this index is valid
1079  if (*p2 != -1)
1080  {
1081  // Check if any of the residuals are larger than the tolerance.
1082  bool all_converged = true;
1083  for (int ii = 0; ii < ensemble_size; ++ii)
1084  {
1085  if (ET::coeff(testvector_[*p2], ii) > ET::coeff(tolerance_, ii))
1086  {
1087  ++ensemble_iterations[ii];
1088  all_converged = false;
1089  }
1090  else if (!(ET::coeff(testvector_[*p2], ii) <= ET::coeff(tolerance_, ii)))
1091  {
1092  // Throw an std::exception if the current residual norm is
1093  // NaN. We know that it's NaN because it is not less than,
1094  // equal to, or greater than the current tolerance. This is
1095  // only possible if either the residual norm or the current
1096  // tolerance is NaN; we assume the former. We also mark the
1097  // test as failed, in case you want to catch the exception.
1098  status_ = Failed;
1099  TEUCHOS_TEST_FOR_EXCEPTION(true, StatusTestError, "StatusTestGenResSubNorm::checkStatus(): NaN has been detected.");
1100  }
1101  }
1102  if (all_converged)
1103  {
1104  ind_[have] = *p2;
1105  have++;
1106  }
1107  }
1108  }
1109  ind_.resize(have);
1110  int need = (quorum_ == -1) ? curNumRHS_ : quorum_;
1111  status_ = (have >= need) ? Passed : Failed;
1112  // Return the current status
1113  return status_;
1114  }
1115 
1117  StatusType getStatus() const { return (status_); };
1119 
1121 
1122 
1124  void reset()
1125  {
1126  status_ = Undefined;
1127  curBlksz_ = 0;
1128  curLSNum_ = 0;
1129  curLSIdx_.resize(0);
1130  numrhs_ = 0;
1131  ind_.resize(0);
1132  firstcallCheckStatus_ = true;
1133  curSoln_ = Teuchos::null;
1134  }
1135 
1137 
1139 
1140 
1142  void print(std::ostream &os, int indent = 0) const
1143  {
1144  os.setf(std::ios_base::scientific);
1145  for (int j = 0; j < indent; j++)
1146  os << ' ';
1147  printStatus(os, status_);
1148  os << resFormStr();
1149  if (status_ == Undefined)
1150  os << ", tol = " << tolerance_ << std::endl;
1151  else
1152  {
1153  os << std::endl;
1154  if (showMaxResNormOnly_ && curBlksz_ > 1)
1155  {
1156  const MagnitudeType maxRelRes = *std::max_element(
1157  testvector_.begin() + curLSIdx_[0], testvector_.begin() + curLSIdx_[curBlksz_ - 1]);
1158  for (int j = 0; j < indent + 13; j++)
1159  os << ' ';
1160  os << "max{residual[" << curLSIdx_[0] << "..." << curLSIdx_[curBlksz_ - 1] << "]} = " << maxRelRes
1161  << (maxRelRes <= tolerance_ ? " <= " : " > ") << tolerance_ << std::endl;
1162  }
1163  else
1164  {
1165  for (int i = 0; i < numrhs_; i++)
1166  {
1167  for (int j = 0; j < indent + 13; j++)
1168  os << ' ';
1169  os << "residual [ " << i << " ] = " << testvector_[i];
1170  os << ((testvector_[i] < tolerance_) ? " < " : (testvector_[i] == tolerance_) ? " == " : (testvector_[i] > tolerance_) ? " > " : " ") << tolerance_ << std::endl;
1171  }
1172  }
1173  }
1174  os << std::endl;
1175  }
1176 
1178  void printStatus(std::ostream &os, StatusType type) const
1179  {
1180  os << std::left << std::setw(13) << std::setfill('.');
1181  switch (type)
1182  {
1183  case Passed:
1184  os << "Converged";
1185  break;
1186  case Failed:
1187  os << "Unconverged";
1188  break;
1189  case Undefined:
1190  default:
1191  os << "**";
1192  break;
1193  }
1194  os << std::left << std::setfill(' ');
1195  return;
1196  }
1198 
1200 
1201 
1203  Teuchos::RCP<MV> getSolution() { return curSoln_; }
1204 
1207  int getQuorum() const { return quorum_; }
1208 
1210  size_t getSubIdx() const { return subIdx_; }
1211 
1213  bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
1214 
1216  std::vector<int> convIndices() { return ind_; }
1217 
1219  MagnitudeType getTolerance() const { return (tolerance_); };
1220 
1222  const std::vector<MagnitudeType> *getTestValue() const { return (&testvector_); };
1223 
1225  const std::vector<MagnitudeType> *getResNormValue() const { return (&resvector_); };
1226 
1228  const std::vector<MagnitudeType> *getScaledNormValue() const { return (&scalevector_); };
1229 
1232  bool getLOADetected() const { return false; }
1233 
1234  const std::vector<int> getEnsembleIterations() const { return ensemble_iterations; }
1235 
1237 
1240 
1246  StatusType firstCallCheckStatusSetup(Iteration<Scalar, MV, OP> *iSolver)
1247  {
1248  int i;
1249  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
1250  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
1251  const LinearProblem<Scalar, MV, OP> &lp = iSolver->getProblem();
1252  // Compute scaling term (done once for each block that's being solved)
1253  if (firstcallCheckStatus_)
1254  {
1255  //
1256  // Get some current solver information.
1257  //
1258  firstcallCheckStatus_ = false;
1259 
1260  // try to access the underlying blocked operator
1261  Teuchos::RCP<const OP> Op = lp.getOperator();
1262  Teuchos::RCP<const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xOp =
1263  Teuchos::rcp_dynamic_cast<const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(Op);
1264  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() << ".");
1265  Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xIntOp =
1266  xOp->getOperator();
1267  TEUCHOS_TEST_FOR_EXCEPTION(xIntOp.is_null(), MueLu::Exceptions::BadCast, "Cannot access Xpetra::Operator stored in Belos::XpetraOperator.");
1268  Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> xMat =
1269  Teuchos::rcp_dynamic_cast<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xIntOp);
1270  TEUCHOS_TEST_FOR_EXCEPTION(xMat.is_null(), MueLu::Exceptions::RuntimeError, "Cannot access Xpetra::Matrix stored in Belos::XpetraOp. Error.");
1271  Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> bMat = Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xMat);
1272  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!");
1273  mapExtractor_ = bMat->getRangeMapExtractor();
1274  TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_.is_null(), MueLu::Exceptions::RuntimeError, "Could not extract map extractor from BlockedCrsMatrix. Error.");
1275  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_ << ".");
1276 
1277  // calculate initial norms
1278  if (scaletype_ == NormOfRHS)
1279  {
1280  Teuchos::RCP<const MV> rhs = lp.getRHS();
1281  numrhs_ = MVT::GetNumberVecs(*rhs);
1282  scalevector_.resize(numrhs_);
1283  MvSubNorm(*rhs, subIdx_, scalevector_, scalenormtype_);
1284  }
1285  else if (scaletype_ == NormOfInitRes)
1286  {
1287  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
1288  numrhs_ = MVT::GetNumberVecs(*init_res);
1289  scalevector_.resize(numrhs_);
1290  MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
1291  }
1292  else if (scaletype_ == NormOfPrecInitRes)
1293  {
1294  Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
1295  numrhs_ = MVT::GetNumberVecs(*init_res);
1296  scalevector_.resize(numrhs_);
1297  MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
1298  }
1299  else if (scaletype_ == NormOfFullInitRes)
1300  {
1301  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
1302  numrhs_ = MVT::GetNumberVecs(*init_res);
1303  scalevector_.resize(numrhs_);
1304  MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
1305  scalevalue_ = one;
1306  }
1307  else if (scaletype_ == NormOfFullPrecInitRes)
1308  {
1309  Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
1310  numrhs_ = MVT::GetNumberVecs(*init_res);
1311  scalevector_.resize(numrhs_);
1312  MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
1313  scalevalue_ = one;
1314  }
1315  else if (scaletype_ == NormOfFullScaledInitRes)
1316  {
1317  Teuchos::RCP<const MV> init_res = lp.getInitResVec();
1318  numrhs_ = MVT::GetNumberVecs(*init_res);
1319  scalevector_.resize(numrhs_);
1320  MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
1321  MvScalingRatio(*init_res, subIdx_, scalevalue_);
1322  }
1323  else if (scaletype_ == NormOfFullScaledPrecInitRes)
1324  {
1325  Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
1326  numrhs_ = MVT::GetNumberVecs(*init_res);
1327  scalevector_.resize(numrhs_);
1328  MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
1329  MvScalingRatio(*init_res, subIdx_, scalevalue_);
1330  }
1331  else
1332  {
1333  numrhs_ = MVT::GetNumberVecs(*(lp.getRHS()));
1334  }
1335 
1336  resvector_.resize(numrhs_);
1337  testvector_.resize(numrhs_);
1338 
1339  curLSNum_ = lp.getLSNumber();
1340  curLSIdx_ = lp.getLSIndex();
1341  curBlksz_ = (int)curLSIdx_.size();
1342  int validLS = 0;
1343  for (i = 0; i < curBlksz_; ++i)
1344  {
1345  if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
1346  validLS++;
1347  }
1348  curNumRHS_ = validLS;
1349  //
1350  // Initialize the testvector.
1351  for (i = 0; i < numrhs_; i++)
1352  {
1353  testvector_[i] = one;
1354  }
1355 
1356  // Return an error if the scaling is zero.
1357  if (scalevalue_ == zero)
1358  {
1359  return Failed;
1360  }
1361  }
1362  return Undefined;
1363  }
1365 
1368 
1370  std::string description() const
1371  {
1372  std::ostringstream oss;
1373  oss << "Belos::StatusTestGenResSubNorm<>: " << resFormStr();
1374  oss << ", tol = " << tolerance_;
1375  return oss.str();
1376  }
1378 
1379 protected:
1380 private:
1382 
1383 
1384  std::string resFormStr() const
1385  {
1386  std::ostringstream oss;
1387  oss << "(";
1388  oss << ((resnormtype_ == OneNorm) ? "1-Norm" : (resnormtype_ == TwoNorm) ? "2-Norm" : "Inf-Norm");
1389  oss << " Exp";
1390  oss << " Res Vec [" << subIdx_ << "]) ";
1391 
1392  // If there is no residual scaling, return current string.
1393  if (scaletype_ != None)
1394  {
1395  // Insert division sign.
1396  oss << "/ ";
1397 
1398  // Determine output string for scaling, if there is any.
1399  if (scaletype_ == UserProvided)
1400  oss << " (User Scale)";
1401  else
1402  {
1403  oss << "(";
1404  oss << ((scalenormtype_ == OneNorm) ? "1-Norm" : (resnormtype_ == TwoNorm) ? "2-Norm" : "Inf-Norm");
1405  if (scaletype_ == NormOfInitRes)
1406  oss << " Res0 [" << subIdx_ << "]";
1407  else if (scaletype_ == NormOfPrecInitRes)
1408  oss << " Prec Res0 [" << subIdx_ << "]";
1409  else if (scaletype_ == NormOfFullInitRes)
1410  oss << " Full Res0 [" << subIdx_ << "]";
1411  else if (scaletype_ == NormOfFullPrecInitRes)
1412  oss << " Full Prec Res0 [" << subIdx_ << "]";
1413  else if (scaletype_ == NormOfFullScaledInitRes)
1414  oss << " scaled Full Res0 [" << subIdx_ << "]";
1415  else if (scaletype_ == NormOfFullScaledPrecInitRes)
1416  oss << " scaled Full Prec Res0 [" << subIdx_ << "]";
1417  else
1418  oss << " RHS [" << subIdx_ << "]";
1419  oss << ")";
1420  }
1421  }
1422 
1423  // TODO add a tagging name
1424 
1425  return oss.str();
1426  }
1427 
1429 
1431 
1432 
1433  // calculate norm of partial multivector
1434  void MvSubNorm(const MV &mv, size_t block, std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normVec, NormType type = TwoNorm)
1435  {
1436 
1437  Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
1438 
1439  Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
1440  typedef MultiVecTraits<Scalar, MV> MVT;
1441  MVT::MvNorm(*SubVec, normVec, type);
1442  }
1443 
1444  // calculate ration of sub-vector length to full vector length (for scalevalue_)
1445  void MvScalingRatio(const MV &mv, size_t block, MagnitudeType &lengthRatio)
1446  {
1447  Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);
1448 
1449  Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
1450 
1451  lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
1452  }
1454 
1456 
1457 
1459  MagnitudeType tolerance_;
1460 
1462  size_t subIdx_;
1463 
1465  int quorum_;
1466 
1468  bool showMaxResNormOnly_;
1469 
1471  NormType resnormtype_;
1472 
1474  ScaleType scaletype_;
1475 
1477  NormType scalenormtype_;
1478 
1480  MagnitudeType scalevalue_;
1481 
1483  std::vector<MagnitudeType> scalevector_;
1484 
1486  std::vector<MagnitudeType> resvector_;
1487 
1489  std::vector<MagnitudeType> testvector_;
1490 
1492  std::vector<int> ind_;
1493 
1495  Teuchos::RCP<MV> curSoln_;
1496 
1498  StatusType status_;
1499 
1501  int curBlksz_;
1502 
1504  int curNumRHS_;
1505 
1507  std::vector<int> curLSIdx_;
1508 
1510  int curLSNum_;
1511 
1513  int numrhs_;
1514 
1516  bool firstcallCheckStatus_;
1517 
1519  bool firstcallDefineResForm_;
1520 
1522  bool firstcallDefineScaleForm_;
1523 
1525  Teuchos::RCP<const ME> mapExtractor_;
1526 
1528  std::vector<int> ensemble_iterations;
1529 
1531 };
1532 
1533 } // namespace Belos
1534 
1535 #endif
1536 
1537 #endif /* BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP */
Belos::StatusTestGenResSubNorm< Scalar, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >, Belos::OperatorT< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > >::MV
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MV
Definition: BelosXpetraStatusTestGenResSubNorm_MP_Vector.hpp:825
Belos::StatusTestGenResSubNorm< Scalar, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >, Belos::OperatorT< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > >::OT
OperatorTraits< Scalar, MV, OP > OT
Definition: BelosXpetraStatusTestGenResSubNorm_MP_Vector.hpp:833
Belos::StatusTestGenResSubNorm< Scalar, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >, Belos::OperatorT< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > >::ME
Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > ME
Definition: BelosXpetraStatusTestGenResSubNorm_MP_Vector.hpp:827
Belos::StatusTestGenResSubNorm< Scalar, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >, Belos::OperatorT< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > >::MagnitudeType
SCT::magnitudeType MagnitudeType
Definition: BelosXpetraStatusTestGenResSubNorm_MP_Vector.hpp:831
EnsembleTraits
Definition: EnsembleTraits.h:8
EnsembleTraits.h
Belos
Definition: BelosStatusTestWeightedGenResNorm_MP_Vector.hpp:80
Belos::StatusTestGenResSubNorm< Scalar, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >, Belos::OperatorT< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > >::MVT
MultiVecTraits< Scalar, MV > MVT
Definition: BelosXpetraStatusTestGenResSubNorm_MP_Vector.hpp:832
Belos::StatusTestGenResSubNorm< Scalar, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >, Belos::OperatorT< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > >::BCRS
Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > BCRS
Definition: BelosXpetraStatusTestGenResSubNorm_MP_Vector.hpp:826
Belos::StatusTestGenResSubNorm< Scalar, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >, Belos::OperatorT< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > >::OP
Belos::OperatorT< MV > OP
Definition: BelosXpetraStatusTestGenResSubNorm_MP_Vector.hpp:828
Belos::StatusTestGenResSubNorm< Scalar, Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >, Belos::OperatorT< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > >::SCT
Teuchos::ScalarTraits< Scalar > SCT
Definition: BelosXpetraStatusTestGenResSubNorm_MP_Vector.hpp:830