Belos  Version of the Day
BelosPseudoBlockCGSolMgr.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
43 #define BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
56 #include "BelosCGSingleRedIter.hpp"
57 #include "BelosCGIter.hpp"
60 #include "BelosStatusTestCombo.hpp"
62 #include "BelosOutputManager.hpp"
63 #include "Teuchos_LAPACK.hpp"
64 #ifdef BELOS_TEUCHOS_TIME_MONITOR
65 #include "Teuchos_TimeMonitor.hpp"
66 #endif
67 
84 namespace Belos {
85 
87 
88 
96  PseudoBlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
97  {}};
98 
106  PseudoBlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
107  {}};
108 
109 
110  // Partial specialization for unsupported ScalarType types.
111  // This contains a stub implementation.
112  template<class ScalarType, class MV, class OP,
113  const bool supportsScalarType =
116  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
117  Belos::Details::LapackSupportsScalar<ScalarType>::value>
118  {
119  static const bool scalarTypeIsSupported =
121  typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
122  scalarTypeIsSupported> base_type;
123 
124  public:
126  base_type ()
127  {}
128  PseudoBlockCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
129  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
130  base_type ()
131  {}
132  virtual ~PseudoBlockCGSolMgr () {}
133 
134  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
135  getResidualStatusTest() const { return Teuchos::null; }
136  };
137 
138 
139  template<class ScalarType, class MV, class OP>
140  class PseudoBlockCGSolMgr<ScalarType, MV, OP, true> :
141  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
142  {
143  private:
146  typedef Teuchos::ScalarTraits<ScalarType> SCT;
147  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
148  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
149 
150  public:
151 
153 
154 
161 
177  PseudoBlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
178  const Teuchos::RCP<Teuchos::ParameterList> &pl );
179 
181  virtual ~PseudoBlockCGSolMgr() {};
182 
184  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
185  return Teuchos::rcp(new PseudoBlockCGSolMgr<ScalarType,MV,OP>);
186  }
188 
190 
191 
192  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
193  return *problem_;
194  }
195 
198  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
199 
202  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
203 
209  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
210  return Teuchos::tuple(timerSolve_);
211  }
212 
213 
224  MagnitudeType achievedTol() const override {
225  return achievedTol_;
226  }
227 
229  int getNumIters() const override {
230  return numIters_;
231  }
232 
236  bool isLOADetected() const override { return false; }
237 
241  ScalarType getConditionEstimate() const {return condEstimate_;}
242 
244  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
245  getResidualStatusTest() const { return convTest_; }
246 
248 
250 
251 
253  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
254 
256  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
257 
259 
261 
262 
266  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
268 
270 
271 
289  ReturnType solve() override;
290 
292 
295 
297  std::string description() const override;
298 
300  private:
301  // Compute the condition number estimate
302  void compute_condnum_tridiag_sym(Teuchos::ArrayView<MagnitudeType> diag,
303  Teuchos::ArrayView<MagnitudeType> offdiag,
304  ScalarType & lambda_min,
305  ScalarType & lambda_max,
306  ScalarType & ConditionNumber );
307 
308  // Linear problem.
309  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
310 
311  // Output manager.
312  Teuchos::RCP<OutputManager<ScalarType> > printer_;
313  Teuchos::RCP<std::ostream> outputStream_;
314 
315  // Status test.
316  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
317  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
318  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
319  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
320 
321  // Current parameter list.
322  Teuchos::RCP<Teuchos::ParameterList> params_;
323 
329  mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
330 
331  // Default solver values.
332  static constexpr int maxIters_default_ = 1000;
333  static constexpr bool assertPositiveDefiniteness_default_ = true;
334  static constexpr bool showMaxResNormOnly_default_ = false;
335  static constexpr int verbosity_default_ = Belos::Errors;
336  static constexpr int outputStyle_default_ = Belos::General;
337  static constexpr int outputFreq_default_ = -1;
338  static constexpr int defQuorum_default_ = 1;
339  static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ = false;
340  static constexpr const char * resScale_default_ = "Norm of Initial Residual";
341  static constexpr const char * label_default_ = "Belos";
342 // https://stackoverflow.com/questions/24398102/constexpr-and-initialization-of-a-static-const-void-pointer-with-reinterpret-cas
343 #if defined(_WIN32) && defined(__clang__)
344  static constexpr std::ostream * outputStream_default_ =
345  __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
346 #else
347  static constexpr std::ostream * outputStream_default_ = &std::cout;
348 #endif
349  static constexpr bool genCondEst_default_ = false;
350 
351  // Current solver values.
352  MagnitudeType convtol_,achievedTol_;
353  int maxIters_, numIters_;
354  int verbosity_, outputStyle_, outputFreq_, defQuorum_;
355  bool assertPositiveDefiniteness_, showMaxResNormOnly_;
356  bool foldConvergenceDetectionIntoAllreduce_;
357  std::string resScale_;
358  bool genCondEst_;
359  ScalarType condEstimate_;
360 
361  // Timers.
362  std::string label_;
363  Teuchos::RCP<Teuchos::Time> timerSolve_;
364 
365  // Internal state variables.
366  bool isSet_;
367  };
368 
369 
370 // Empty Constructor
371 template<class ScalarType, class MV, class OP>
373  outputStream_(Teuchos::rcp(outputStream_default_,false)),
374  convtol_(DefaultSolverParameters::convTol),
375  maxIters_(maxIters_default_),
376  numIters_(0),
377  verbosity_(verbosity_default_),
378  outputStyle_(outputStyle_default_),
379  outputFreq_(outputFreq_default_),
380  defQuorum_(defQuorum_default_),
381  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
382  showMaxResNormOnly_(showMaxResNormOnly_default_),
383  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
384  resScale_(resScale_default_),
385  genCondEst_(genCondEst_default_),
386  condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
387  label_(label_default_),
388  isSet_(false)
389 {}
390 
391 // Basic Constructor
392 template<class ScalarType, class MV, class OP>
395  const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
396  problem_(problem),
397  outputStream_(Teuchos::rcp(outputStream_default_,false)),
398  convtol_(DefaultSolverParameters::convTol),
399  maxIters_(maxIters_default_),
400  numIters_(0),
401  verbosity_(verbosity_default_),
402  outputStyle_(outputStyle_default_),
403  outputFreq_(outputFreq_default_),
404  defQuorum_(defQuorum_default_),
405  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
406  showMaxResNormOnly_(showMaxResNormOnly_default_),
407  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
408  resScale_(resScale_default_),
409  genCondEst_(genCondEst_default_),
410  condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
411  label_(label_default_),
412  isSet_(false)
413 {
414  TEUCHOS_TEST_FOR_EXCEPTION(
415  problem_.is_null (), std::invalid_argument,
416  "Belos::PseudoBlockCGSolMgr two-argument constructor: "
417  "'problem' is null. You must supply a non-null Belos::LinearProblem "
418  "instance when calling this constructor.");
419 
420  if (! pl.is_null ()) {
421  // Set the parameters using the list that was passed in.
422  setParameters (pl);
423  }
424 }
425 
426 template<class ScalarType, class MV, class OP>
427 void
429 setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
430 {
431  using Teuchos::ParameterList;
432  using Teuchos::parameterList;
433  using Teuchos::RCP;
434  using Teuchos::rcp;
435 
436  RCP<const ParameterList> defaultParams = this->getValidParameters ();
437 
438  // Create the internal parameter list if one doesn't already exist.
439  // Belos' solvers treat the input ParameterList to setParameters as
440  // a "delta" -- that is, a change from the current state -- so the
441  // default parameter list (if the input is null) should be empty.
442  // This explains also why Belos' solvers copy parameters one by one
443  // from the input list to the current list.
444  //
445  // Belos obfuscates the latter, because it takes the input parameter
446  // list by RCP, rather than by (nonconst) reference. The latter
447  // would make more sense, given that it doesn't actually keep the
448  // input parameter list.
449  //
450  // Note, however, that Belos still correctly triggers the "used"
451  // field of each parameter in the input list. While isParameter()
452  // doesn't (apparently) trigger the "used" flag, get() certainly
453  // does.
454 
455  if (params_.is_null ()) {
456  // Create an empty list with the same name as the default list.
457  params_ = parameterList (defaultParams->name ());
458  } else {
459  params->validateParameters (*defaultParams);
460  }
461 
462  // Check for maximum number of iterations
463  if (params->isParameter ("Maximum Iterations")) {
464  maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
465 
466  // Update parameter in our list and in status test.
467  params_->set ("Maximum Iterations", maxIters_);
468  if (! maxIterTest_.is_null ()) {
469  maxIterTest_->setMaxIters (maxIters_);
470  }
471  }
472 
473  // Check if positive definiteness assertions are to be performed
474  if (params->isParameter ("Assert Positive Definiteness")) {
475  assertPositiveDefiniteness_ =
476  params->get ("Assert Positive Definiteness",
477  assertPositiveDefiniteness_default_);
478 
479  // Update parameter in our list.
480  params_->set ("Assert Positive Definiteness", assertPositiveDefiniteness_);
481  }
482 
483  if (params->isParameter("Fold Convergence Detection Into Allreduce")) {
484  foldConvergenceDetectionIntoAllreduce_ = params->get("Fold Convergence Detection Into Allreduce",
485  foldConvergenceDetectionIntoAllreduce_default_);
486  }
487 
488  // Check to see if the timer label changed.
489  if (params->isParameter ("Timer Label")) {
490  const std::string tempLabel = params->get ("Timer Label", label_default_);
491 
492  // Update parameter in our list and solver timer
493  if (tempLabel != label_) {
494  label_ = tempLabel;
495  params_->set ("Timer Label", label_);
496  const std::string solveLabel =
497  label_ + ": PseudoBlockCGSolMgr total solve time";
498 #ifdef BELOS_TEUCHOS_TIME_MONITOR
499  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
500 #endif
501  }
502  }
503 
504  // Check for a change in verbosity level
505  if (params->isParameter ("Verbosity")) {
506  if (Teuchos::isParameterType<int> (*params, "Verbosity")) {
507  verbosity_ = params->get ("Verbosity", verbosity_default_);
508  } else {
509  verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params, "Verbosity");
510  }
511 
512  // Update parameter in our list.
513  params_->set ("Verbosity", verbosity_);
514  if (! printer_.is_null ()) {
515  printer_->setVerbosity (verbosity_);
516  }
517  }
518 
519  // Check for a change in output style
520  if (params->isParameter ("Output Style")) {
521  if (Teuchos::isParameterType<int> (*params, "Output Style")) {
522  outputStyle_ = params->get ("Output Style", outputStyle_default_);
523  } else {
524  // FIXME (mfh 29 Jul 2015) What if the type is wrong?
525  outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params, "Output Style");
526  }
527 
528  // Reconstruct the convergence test if the explicit residual test
529  // is not being used.
530  params_->set ("Output Style", outputStyle_);
531  outputTest_ = Teuchos::null;
532  }
533 
534  // output stream
535  if (params->isParameter ("Output Stream")) {
536  outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
537 
538  // Update parameter in our list.
539  params_->set ("Output Stream", outputStream_);
540  if (! printer_.is_null ()) {
541  printer_->setOStream (outputStream_);
542  }
543  }
544 
545  // frequency level
546  if (verbosity_ & Belos::StatusTestDetails) {
547  if (params->isParameter ("Output Frequency")) {
548  outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
549  }
550 
551  // Update parameter in out list and output status test.
552  params_->set ("Output Frequency", outputFreq_);
553  if (! outputTest_.is_null ()) {
554  outputTest_->setOutputFrequency (outputFreq_);
555  }
556  }
557 
558  // Condition estimate
559  if (params->isParameter ("Estimate Condition Number")) {
560  genCondEst_ = params->get ("Estimate Condition Number", genCondEst_default_);
561  }
562 
563  // Create output manager if we need to.
564  if (printer_.is_null ()) {
565  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
566  }
567 
568  // Convergence
569  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
570  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
571 
572  // Check for convergence tolerance
573  if (params->isParameter ("Convergence Tolerance")) {
574  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
575  convtol_ = params->get ("Convergence Tolerance",
576  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
577  }
578  else {
579  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
580  }
581 
582  // Update parameter in our list and residual tests.
583  params_->set ("Convergence Tolerance", convtol_);
584  if (! convTest_.is_null ()) {
585  convTest_->setTolerance (convtol_);
586  }
587  }
588 
589  if (params->isParameter ("Show Maximum Residual Norm Only")) {
590  showMaxResNormOnly_ = params->get<bool> ("Show Maximum Residual Norm Only");
591 
592  // Update parameter in our list and residual tests
593  params_->set ("Show Maximum Residual Norm Only", showMaxResNormOnly_);
594  if (! convTest_.is_null ()) {
595  convTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
596  }
597  }
598 
599  // Check for a change in scaling, if so we need to build new residual tests.
600  bool newResTest = false;
601  {
602  // "Residual Scaling" is the old parameter name; "Implicit
603  // Residual Scaling" is the new name. We support both options for
604  // backwards compatibility.
605  std::string tempResScale = resScale_;
606  bool implicitResidualScalingName = false;
607  if (params->isParameter ("Residual Scaling")) {
608  tempResScale = params->get<std::string> ("Residual Scaling");
609  }
610  else if (params->isParameter ("Implicit Residual Scaling")) {
611  tempResScale = params->get<std::string> ("Implicit Residual Scaling");
612  implicitResidualScalingName = true;
613  }
614 
615  // Only update the scaling if it's different.
616  if (resScale_ != tempResScale) {
617  const Belos::ScaleType resScaleType =
618  convertStringToScaleType (tempResScale);
619  resScale_ = tempResScale;
620 
621  // Update parameter in our list and residual tests, using the
622  // given parameter name.
623  if (implicitResidualScalingName) {
624  params_->set ("Implicit Residual Scaling", resScale_);
625  }
626  else {
627  params_->set ("Residual Scaling", resScale_);
628  }
629 
630  if (! convTest_.is_null ()) {
631  try {
632  convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
633  }
634  catch (std::exception& e) {
635  // Make sure the convergence test gets constructed again.
636  newResTest = true;
637  }
638  }
639  }
640  }
641 
642  // Get the deflation quorum, or number of converged systems before deflation is allowed
643  if (params->isParameter ("Deflation Quorum")) {
644  defQuorum_ = params->get ("Deflation Quorum", defQuorum_);
645  params_->set ("Deflation Quorum", defQuorum_);
646  if (! convTest_.is_null ()) {
647  convTest_->setQuorum( defQuorum_ );
648  }
649  }
650 
651  // Create status tests if we need to.
652 
653  // Basic test checks maximum iterations and native residual.
654  if (maxIterTest_.is_null ()) {
655  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
656  }
657 
658  // Implicit residual test, using the native residual to determine if convergence was achieved.
659  if (convTest_.is_null () || newResTest) {
660  convTest_ = rcp (new StatusTestResNorm_t (convtol_, defQuorum_, showMaxResNormOnly_));
661  convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
662  }
663 
664  if (sTest_.is_null () || newResTest) {
665  sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
666  }
667 
668  if (outputTest_.is_null () || newResTest) {
669  // Create the status test output class.
670  // This class manages and formats the output from the status test.
671  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
672  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
674 
675  // Set the solver string for the output test
676  const std::string solverDesc = " Pseudo Block CG ";
677  outputTest_->setSolverDesc (solverDesc);
678  }
679 
680  // Create the timer if we need to.
681  if (timerSolve_.is_null ()) {
682  const std::string solveLabel =
683  label_ + ": PseudoBlockCGSolMgr total solve time";
684 #ifdef BELOS_TEUCHOS_TIME_MONITOR
685  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
686 #endif
687  }
688 
689  // Inform the solver manager that the current parameters were set.
690  isSet_ = true;
691 }
692 
693 
694 template<class ScalarType, class MV, class OP>
695 Teuchos::RCP<const Teuchos::ParameterList>
697 {
698  using Teuchos::ParameterList;
699  using Teuchos::parameterList;
700  using Teuchos::RCP;
701 
702  if (validParams_.is_null()) {
703  // Set all the valid parameters and their default values.
704  RCP<ParameterList> pl = parameterList ();
705  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
706  "The relative residual tolerance that needs to be achieved by the\n"
707  "iterative solver in order for the linear system to be declared converged.");
708  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
709  "The maximum number of block iterations allowed for each\n"
710  "set of RHS solved.");
711  pl->set("Assert Positive Definiteness", static_cast<bool>(assertPositiveDefiniteness_default_),
712  "Whether or not to assert that the linear operator\n"
713  "and the preconditioner are indeed positive definite.");
714  pl->set("Verbosity", static_cast<int>(verbosity_default_),
715  "What type(s) of solver information should be outputted\n"
716  "to the output stream.");
717  pl->set("Output Style", static_cast<int>(outputStyle_default_),
718  "What style is used for the solver information outputted\n"
719  "to the output stream.");
720  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
721  "How often convergence information should be outputted\n"
722  "to the output stream.");
723  pl->set("Deflation Quorum", static_cast<int>(defQuorum_default_),
724  "The number of linear systems that need to converge before\n"
725  "they are deflated. This number should be <= block size.");
726  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
727  "A reference-counted pointer to the output stream where all\n"
728  "solver output is sent.");
729  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
730  "When convergence information is printed, only show the maximum\n"
731  "relative residual norm when the block size is greater than one.");
732  pl->set("Implicit Residual Scaling", resScale_default_,
733  "The type of scaling used in the residual convergence test.");
734  pl->set("Estimate Condition Number", static_cast<bool>(genCondEst_default_),
735  "Whether or not to estimate the condition number of the preconditioned system.");
736  // We leave the old name as a valid parameter for backwards
737  // compatibility (so that validateParametersAndSetDefaults()
738  // doesn't raise an exception if it encounters "Residual
739  // Scaling"). The new name was added for compatibility with other
740  // solvers, none of which use "Residual Scaling".
741  pl->set("Residual Scaling", resScale_default_,
742  "The type of scaling used in the residual convergence test. This "
743  "name is deprecated; the new name is \"Implicit Residual Scaling\".");
744  pl->set("Timer Label", static_cast<const char *>(label_default_),
745  "The string to use as a prefix for the timer labels.");
746  pl->set("Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
747  "Merge the allreduce for convergence detection with the one for CG.\n"
748  "This saves one all-reduce, but incurs more computation.");
749  validParams_ = pl;
750  }
751  return validParams_;
752 }
753 
754 
755 // solve()
756 template<class ScalarType, class MV, class OP>
758 {
759  const char prefix[] = "Belos::PseudoBlockCGSolMgr::solve: ";
760 
761  // Set the current parameters if they were not set before.
762  // NOTE: This may occur if the user generated the solver manager with the default constructor and
763  // then didn't set any parameters using setParameters().
764  if (!isSet_) { setParameters( params_ ); }
765 
766  TEUCHOS_TEST_FOR_EXCEPTION
767  (! problem_->isProblemSet (), PseudoBlockCGSolMgrLinearProblemFailure,
768  prefix << "The linear problem to solve is not ready. You must call "
769  "setProblem() on the Belos::LinearProblem instance before telling the "
770  "Belos solver to solve it.");
771 
772  // Create indices for the linear systems to be solved.
773  int startPtr = 0;
774  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
775  int numCurrRHS = numRHS2Solve;
776 
777  std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
778  for (int i=0; i<numRHS2Solve; ++i) {
779  currIdx[i] = startPtr+i;
780  currIdx2[i]=i;
781  }
782 
783  // Inform the linear problem of the current linear system to solve.
784  problem_->setLSIndex( currIdx );
785 
787  // Parameter list (iteration)
788  Teuchos::ParameterList plist;
789 
790  plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
791  if(genCondEst_) plist.set("Max Size For Condest",maxIters_);
792 
793  // Reset the status test.
794  outputTest_->reset();
795 
796  // Assume convergence is achieved, then let any failed convergence set this to false.
797  bool isConverged = true;
798 
800  // Pseudo-Block CG solver
801  Teuchos::RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
802  if (numRHS2Solve == 1) {
803  plist.set("Fold Convergence Detection Into Allreduce",
804  foldConvergenceDetectionIntoAllreduce_);
805  block_cg_iter =
806  Teuchos::rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, convTest_, plist));
807  } else {
808  block_cg_iter =
809  Teuchos::rcp (new PseudoBlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
810  }
811 
812  // Setup condition estimate
813  block_cg_iter->setDoCondEst(genCondEst_);
814  bool condEstPerf = false;
815 
816  // Enter solve() iterations
817  {
818 #ifdef BELOS_TEUCHOS_TIME_MONITOR
819  Teuchos::TimeMonitor slvtimer(*timerSolve_);
820 #endif
821 
822  while ( numRHS2Solve > 0 ) {
823 
824  // Reset the active / converged vectors from this block
825  std::vector<int> convRHSIdx;
826  std::vector<int> currRHSIdx( currIdx );
827  currRHSIdx.resize(numCurrRHS);
828 
829  // Reset the number of iterations.
830  block_cg_iter->resetNumIters();
831 
832  // Reset the number of calls that the status test output knows about.
833  outputTest_->resetNumCalls();
834 
835  // Get the current residual for this block of linear systems.
836  Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
837 
838  // Get a new state struct and initialize the solver.
840  newState.R = R_0;
841  block_cg_iter->initializeCG(newState);
842 
843  while(1) {
844 
845  // tell block_gmres_iter to iterate
846  try {
847 
848  block_cg_iter->iterate();
849 
851  //
852  // check convergence first
853  //
855  if ( convTest_->getStatus() == Passed ) {
856 
857  // Figure out which linear systems converged.
858  std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
859 
860  // If the number of converged linear systems is equal to the
861  // number of current linear systems, then we are done with this block.
862  if (convIdx.size() == currRHSIdx.size())
863  break; // break from while(1){block_cg_iter->iterate()}
864 
865  // Inform the linear problem that we are finished with this current linear system.
866  problem_->setCurrLS();
867 
868  // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
869  int have = 0;
870  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
871  bool found = false;
872  for (unsigned int j=0; j<convIdx.size(); ++j) {
873  if (currRHSIdx[i] == convIdx[j]) {
874  found = true;
875  break;
876  }
877  }
878  if (!found) {
879  currIdx2[have] = currIdx2[i];
880  currRHSIdx[have++] = currRHSIdx[i];
881  }
882  }
883  currRHSIdx.resize(have);
884  currIdx2.resize(have);
885 
886  // Compute condition estimate if the very first linear system in the block has converged.
887  if (currRHSIdx[0] != 0 && genCondEst_ && !condEstPerf)
888  {
889  // Compute the estimate.
890  ScalarType l_min, l_max;
891  Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
892  Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
893  compute_condnum_tridiag_sym(diag,offdiag,l_min,l_max,condEstimate_);
894 
895  // Make sure not to do more condition estimate computations for this solve.
896  block_cg_iter->setDoCondEst(false);
897  condEstPerf = true;
898  }
899 
900  // Set the remaining indices after deflation.
901  problem_->setLSIndex( currRHSIdx );
902 
903  // Get the current residual vector.
904  std::vector<MagnitudeType> norms;
905  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
906  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
907 
908  // Set the new state and initialize the solver.
910  defstate.R = R_0;
911  block_cg_iter->initializeCG(defstate);
912  }
913 
915  //
916  // check for maximum iterations
917  //
919  else if ( maxIterTest_->getStatus() == Passed ) {
920  // we don't have convergence
921  isConverged = false;
922  break; // break from while(1){block_cg_iter->iterate()}
923  }
924 
926  //
927  // we returned from iterate(), but none of our status tests Passed.
928  // something is wrong, and it is probably our fault.
929  //
931 
932  else {
933  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
934  "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
935  }
936  }
937  catch (const std::exception &e) {
938  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
939  << block_cg_iter->getNumIters() << std::endl
940  << e.what() << std::endl;
941  throw;
942  }
943  }
944 
945  // Inform the linear problem that we are finished with this block linear system.
946  problem_->setCurrLS();
947 
948  // Update indices for the linear systems to be solved.
949  startPtr += numCurrRHS;
950  numRHS2Solve -= numCurrRHS;
951 
952  if ( numRHS2Solve > 0 ) {
953 
954  numCurrRHS = numRHS2Solve;
955  currIdx.resize( numCurrRHS );
956  currIdx2.resize( numCurrRHS );
957  for (int i=0; i<numCurrRHS; ++i)
958  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
959 
960  // Set the next indices.
961  problem_->setLSIndex( currIdx );
962  }
963  else {
964  currIdx.resize( numRHS2Solve );
965  }
966 
967  }// while ( numRHS2Solve > 0 )
968 
969  }
970 
971  // print final summary
972  sTest_->print( printer_->stream(FinalSummary) );
973 
974  // print timing information
975 #ifdef BELOS_TEUCHOS_TIME_MONITOR
976  // Calling summarize() can be expensive, so don't call unless the
977  // user wants to print out timing details. summarize() will do all
978  // the work even if it's passed a "black hole" output stream.
979  if (verbosity_ & TimingDetails)
980  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
981 #endif
982 
983  // get iteration information for this solve
984  numIters_ = maxIterTest_->getNumIters();
985 
986  // Save the convergence test value ("achieved tolerance") for this
987  // solve.
988  const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
989  if (pTestValues != NULL && pTestValues->size () > 0) {
990  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
991  }
992 
993  // Do condition estimate, if needed
994  if (genCondEst_ && !condEstPerf) {
995  ScalarType l_min, l_max;
996  Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
997  Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
998  compute_condnum_tridiag_sym(diag,offdiag,l_min,l_max,condEstimate_);
999  condEstPerf = true;
1000  }
1001 
1002  if (! isConverged) {
1003  return Unconverged; // return from PseudoBlockCGSolMgr::solve()
1004  }
1005  return Converged; // return from PseudoBlockCGSolMgr::solve()
1006 }
1007 
1008 // This method requires the solver manager to return a std::string that describes itself.
1009 template<class ScalarType, class MV, class OP>
1011 {
1012  std::ostringstream oss;
1013  oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1014  oss << "{";
1015  oss << "}";
1016  return oss.str();
1017 }
1018 
1019 
1020 template<class ScalarType, class MV, class OP>
1021 void
1023 compute_condnum_tridiag_sym (Teuchos::ArrayView<MagnitudeType> diag,
1024  Teuchos::ArrayView<MagnitudeType> offdiag,
1025  ScalarType & lambda_min,
1026  ScalarType & lambda_max,
1027  ScalarType & ConditionNumber )
1028 {
1029  typedef Teuchos::ScalarTraits<ScalarType> STS;
1030 
1031  /* Copied from az_cg.c: compute_condnum_tridiag_sym */
1032  /* diag == ScalarType vector of size N, containing the diagonal
1033  elements of A
1034  offdiag == ScalarType vector of size N-1, containing the offdiagonal
1035  elements of A. Note that A is supposed to be symmatric
1036  */
1037  int info = 0;
1038  const int N = diag.size ();
1039  ScalarType scalar_dummy;
1040  std::vector<MagnitudeType> mag_dummy(4*N);
1041  char char_N = 'N';
1042  Teuchos::LAPACK<int,ScalarType> lapack;
1043 
1044  lambda_min = STS::one ();
1045  lambda_max = STS::one ();
1046  if( N > 2 ) {
1047  lapack.PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
1048  &scalar_dummy, 1, &mag_dummy[0], &info);
1049  TEUCHOS_TEST_FOR_EXCEPTION
1050  (info < 0, std::logic_error, "Belos::PseudoBlockCGSolMgr::"
1051  "compute_condnum_tridiag_sym: LAPACK's _PTEQR failed with info = "
1052  << info << " < 0. This suggests there might be a bug in the way Belos "
1053  "is calling LAPACK. Please report this to the Belos developers.");
1054  lambda_min = Teuchos::as<ScalarType> (diag[N-1]);
1055  lambda_max = Teuchos::as<ScalarType> (diag[0]);
1056  }
1057 
1058  // info > 0 means that LAPACK's eigensolver didn't converge. This
1059  // is unlikely but may be possible. In that case, the best we can
1060  // do is use the eigenvalues that it computes, as long as lambda_max
1061  // >= lambda_min.
1062  if (STS::real (lambda_max) < STS::real (lambda_min)) {
1063  ConditionNumber = STS::one ();
1064  }
1065  else {
1066  // It's OK for the condition number to be Inf.
1067  ConditionNumber = lambda_max / lambda_min;
1068  }
1069 
1070 } /* compute_condnum_tridiag_sym */
1071 
1072 
1073 
1074 
1075 
1076 } // end Belos namespace
1077 
1078 #endif /* BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP */
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
ScalarType getConditionEstimate() const
Gets the estimated condition number.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:120
PseudoBlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:79
Structure to contain pointers to CGIteration state variables.
Belos concrete class for performing the conjugate-gradient (CG) iteration.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
A factory class for generating StatusTestOutput objects.
PseudoBlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
An implementation of StatusTestResNorm using a family of residual norms.
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:293
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
A factory class for generating StatusTestOutput objects.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Traits class which defines basic operations on multivectors.
The Belos::PseudoBlockCGSolMgr provides a powerful and fully-featured solver manager over the pseudo-...
Belos::StatusTest for logically combining several status tests.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
Pure virtual base class which describes the basic interface for a solver manager. ...
Belos concrete class for performing a single-reduction conjugate-gradient (CG) iteration.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
PseudoBlockCGSolMgrOrthoFailure(const std::string &what_arg)
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
Return the residual status test.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
PseudoBlockCGSolMgrLinearProblemFailure(const std::string &what_arg)
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
A class for extending the status testing capabilities of Belos via logical combinations.
PseudoBlockCGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate or...
Class which defines basic traits for the operator type.
Belos concrete class for performing the pseudo-block CG iteration.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:283
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.

Generated for Belos by doxygen 1.8.14