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_BLAS.hpp"
64 #include "Teuchos_LAPACK.hpp"
65 #ifdef BELOS_TEUCHOS_TIME_MONITOR
66 #include "Teuchos_TimeMonitor.hpp"
67 #endif
68 
85 namespace Belos {
86 
88 
89 
97  PseudoBlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
98  {}};
99 
107  PseudoBlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
108  {}};
109 
110 
111  // Partial specialization for unsupported ScalarType types.
112  // This contains a stub implementation.
113  template<class ScalarType, class MV, class OP,
114  const bool supportsScalarType =
117  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
118  Belos::Details::LapackSupportsScalar<ScalarType>::value>
119  {
120  static const bool scalarTypeIsSupported =
122  typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
123  scalarTypeIsSupported> base_type;
124 
125  public:
127  base_type ()
128  {}
129  PseudoBlockCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
130  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
131  base_type ()
132  {}
133  virtual ~PseudoBlockCGSolMgr () {}
134 
135  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
136  getResidualStatusTest() const { return Teuchos::null; }
137  };
138 
139 
140  template<class ScalarType, class MV, class OP>
141  class PseudoBlockCGSolMgr<ScalarType, MV, OP, true> :
142  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
143  {
144  private:
147  typedef Teuchos::ScalarTraits<ScalarType> SCT;
148  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
149  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
150 
151  public:
152 
154 
155 
162 
178  PseudoBlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
179  const Teuchos::RCP<Teuchos::ParameterList> &pl );
180 
182  virtual ~PseudoBlockCGSolMgr() {};
183 
185  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
186  return Teuchos::rcp(new PseudoBlockCGSolMgr<ScalarType,MV,OP>);
187  }
189 
191 
192 
193  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
194  return *problem_;
195  }
196 
199  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
200 
203  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
204 
210  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
211  return Teuchos::tuple(timerSolve_);
212  }
213 
214 
225  MagnitudeType achievedTol() const override {
226  return achievedTol_;
227  }
228 
230  int getNumIters() const override {
231  return numIters_;
232  }
233 
237  bool isLOADetected() const override { return false; }
238 
242  ScalarType getConditionEstimate() const {return condEstimate_;}
243 
245  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> >
246  getResidualStatusTest() const { return convTest_; }
247 
249 
251 
252 
254  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
255 
257  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
258 
260 
262 
263 
267  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
269 
271 
272 
290  ReturnType solve() override;
291 
293 
296 
298  std::string description() const override;
299 
301  private:
302  // Compute the condition number estimate
303  void compute_condnum_tridiag_sym(Teuchos::ArrayView<MagnitudeType> diag,
304  Teuchos::ArrayView<MagnitudeType> offdiag,
305  ScalarType & lambda_min,
306  ScalarType & lambda_max,
307  ScalarType & ConditionNumber );
308 
309  // Linear problem.
310  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
311 
312  // Output manager.
313  Teuchos::RCP<OutputManager<ScalarType> > printer_;
314  Teuchos::RCP<std::ostream> outputStream_;
315 
316  // Status test.
317  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
318  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
319  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
320  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
321 
322  // Current parameter list.
323  Teuchos::RCP<Teuchos::ParameterList> params_;
324 
330  mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
331 
332  // Default solver values.
333  static constexpr int maxIters_default_ = 1000;
334  static constexpr bool assertPositiveDefiniteness_default_ = true;
335  static constexpr bool showMaxResNormOnly_default_ = false;
336  static constexpr int verbosity_default_ = Belos::Errors;
337  static constexpr int outputStyle_default_ = Belos::General;
338  static constexpr int outputFreq_default_ = -1;
339  static constexpr int defQuorum_default_ = 1;
340  static constexpr const char * resScale_default_ = "Norm of Initial Residual";
341  static constexpr const char * label_default_ = "Belos";
342  static constexpr std::ostream * outputStream_default_ = &std::cout;
343  static constexpr bool genCondEst_default_ = false;
344 
345  // Current solver values.
346  MagnitudeType convtol_,achievedTol_;
347  int maxIters_, numIters_;
348  int verbosity_, outputStyle_, outputFreq_, defQuorum_;
349  bool assertPositiveDefiniteness_, showMaxResNormOnly_;
350  std::string resScale_;
351  bool genCondEst_;
352  ScalarType condEstimate_;
353 
354  // Timers.
355  std::string label_;
356  Teuchos::RCP<Teuchos::Time> timerSolve_;
357 
358  // Internal state variables.
359  bool isSet_;
360  };
361 
362 
363 // Empty Constructor
364 template<class ScalarType, class MV, class OP>
366  outputStream_(Teuchos::rcp(outputStream_default_,false)),
367  convtol_(DefaultSolverParameters::convTol),
368  maxIters_(maxIters_default_),
369  numIters_(0),
370  verbosity_(verbosity_default_),
371  outputStyle_(outputStyle_default_),
372  outputFreq_(outputFreq_default_),
373  defQuorum_(defQuorum_default_),
374  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
375  showMaxResNormOnly_(showMaxResNormOnly_default_),
376  resScale_(resScale_default_),
377  genCondEst_(genCondEst_default_),
378  condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
379  label_(label_default_),
380  isSet_(false)
381 {}
382 
383 // Basic Constructor
384 template<class ScalarType, class MV, class OP>
387  const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
388  problem_(problem),
389  outputStream_(Teuchos::rcp(outputStream_default_,false)),
390  convtol_(DefaultSolverParameters::convTol),
391  maxIters_(maxIters_default_),
392  numIters_(0),
393  verbosity_(verbosity_default_),
394  outputStyle_(outputStyle_default_),
395  outputFreq_(outputFreq_default_),
396  defQuorum_(defQuorum_default_),
397  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
398  showMaxResNormOnly_(showMaxResNormOnly_default_),
399  resScale_(resScale_default_),
400  genCondEst_(genCondEst_default_),
401  condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
402  label_(label_default_),
403  isSet_(false)
404 {
405  TEUCHOS_TEST_FOR_EXCEPTION(
406  problem_.is_null (), std::invalid_argument,
407  "Belos::PseudoBlockCGSolMgr two-argument constructor: "
408  "'problem' is null. You must supply a non-null Belos::LinearProblem "
409  "instance when calling this constructor.");
410 
411  if (! pl.is_null ()) {
412  // Set the parameters using the list that was passed in.
413  setParameters (pl);
414  }
415 }
416 
417 template<class ScalarType, class MV, class OP>
418 void
420 setParameters (const Teuchos::RCP<Teuchos::ParameterList>& params)
421 {
422  using Teuchos::ParameterList;
423  using Teuchos::parameterList;
424  using Teuchos::RCP;
425  using Teuchos::rcp;
426 
427  RCP<const ParameterList> defaultParams = this->getValidParameters ();
428 
429  // Create the internal parameter list if one doesn't already exist.
430  // Belos' solvers treat the input ParameterList to setParameters as
431  // a "delta" -- that is, a change from the current state -- so the
432  // default parameter list (if the input is null) should be empty.
433  // This explains also why Belos' solvers copy parameters one by one
434  // from the input list to the current list.
435  //
436  // Belos obfuscates the latter, because it takes the input parameter
437  // list by RCP, rather than by (nonconst) reference. The latter
438  // would make more sense, given that it doesn't actually keep the
439  // input parameter list.
440  //
441  // Note, however, that Belos still correctly triggers the "used"
442  // field of each parameter in the input list. While isParameter()
443  // doesn't (apparently) trigger the "used" flag, get() certainly
444  // does.
445 
446  if (params_.is_null ()) {
447  // Create an empty list with the same name as the default list.
448  params_ = parameterList (defaultParams->name ());
449  } else {
450  params->validateParameters (*defaultParams);
451  }
452 
453  // Check for maximum number of iterations
454  if (params->isParameter ("Maximum Iterations")) {
455  maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
456 
457  // Update parameter in our list and in status test.
458  params_->set ("Maximum Iterations", maxIters_);
459  if (! maxIterTest_.is_null ()) {
460  maxIterTest_->setMaxIters (maxIters_);
461  }
462  }
463 
464  // Check if positive definiteness assertions are to be performed
465  if (params->isParameter ("Assert Positive Definiteness")) {
466  assertPositiveDefiniteness_ =
467  params->get ("Assert Positive Definiteness",
468  assertPositiveDefiniteness_default_);
469 
470  // Update parameter in our list.
471  params_->set ("Assert Positive Definiteness", assertPositiveDefiniteness_);
472  }
473 
474  // Check to see if the timer label changed.
475  if (params->isParameter ("Timer Label")) {
476  const std::string tempLabel = params->get ("Timer Label", label_default_);
477 
478  // Update parameter in our list and solver timer
479  if (tempLabel != label_) {
480  label_ = tempLabel;
481  params_->set ("Timer Label", label_);
482  const std::string solveLabel =
483  label_ + ": PseudoBlockCGSolMgr total solve time";
484 #ifdef BELOS_TEUCHOS_TIME_MONITOR
485  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
486 #endif
487  }
488  }
489 
490  // Check for a change in verbosity level
491  if (params->isParameter ("Verbosity")) {
492  if (Teuchos::isParameterType<int> (*params, "Verbosity")) {
493  verbosity_ = params->get ("Verbosity", verbosity_default_);
494  } else {
495  verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params, "Verbosity");
496  }
497 
498  // Update parameter in our list.
499  params_->set ("Verbosity", verbosity_);
500  if (! printer_.is_null ()) {
501  printer_->setVerbosity (verbosity_);
502  }
503  }
504 
505  // Check for a change in output style
506  if (params->isParameter ("Output Style")) {
507  if (Teuchos::isParameterType<int> (*params, "Output Style")) {
508  outputStyle_ = params->get ("Output Style", outputStyle_default_);
509  } else {
510  // FIXME (mfh 29 Jul 2015) What if the type is wrong?
511  outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params, "Output Style");
512  }
513 
514  // Reconstruct the convergence test if the explicit residual test
515  // is not being used.
516  params_->set ("Output Style", outputStyle_);
517  outputTest_ = Teuchos::null;
518  }
519 
520  // output stream
521  if (params->isParameter ("Output Stream")) {
522  outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
523 
524  // Update parameter in our list.
525  params_->set ("Output Stream", outputStream_);
526  if (! printer_.is_null ()) {
527  printer_->setOStream (outputStream_);
528  }
529  }
530 
531  // frequency level
532  if (verbosity_ & Belos::StatusTestDetails) {
533  if (params->isParameter ("Output Frequency")) {
534  outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
535  }
536 
537  // Update parameter in out list and output status test.
538  params_->set ("Output Frequency", outputFreq_);
539  if (! outputTest_.is_null ()) {
540  outputTest_->setOutputFrequency (outputFreq_);
541  }
542  }
543 
544  // Condition estimate
545  if (params->isParameter ("Estimate Condition Number")) {
546  genCondEst_ = params->get ("Estimate Condition Number", genCondEst_default_);
547  }
548 
549  // Create output manager if we need to.
550  if (printer_.is_null ()) {
551  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
552  }
553 
554  // Convergence
555  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
556  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
557 
558  // Check for convergence tolerance
559  if (params->isParameter ("Convergence Tolerance")) {
560  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
561  convtol_ = params->get ("Convergence Tolerance",
562  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
563  }
564  else {
565  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
566  }
567 
568  // Update parameter in our list and residual tests.
569  params_->set ("Convergence Tolerance", convtol_);
570  if (! convTest_.is_null ()) {
571  convTest_->setTolerance (convtol_);
572  }
573  }
574 
575  if (params->isParameter ("Show Maximum Residual Norm Only")) {
576  showMaxResNormOnly_ = params->get<bool> ("Show Maximum Residual Norm Only");
577 
578  // Update parameter in our list and residual tests
579  params_->set ("Show Maximum Residual Norm Only", showMaxResNormOnly_);
580  if (! convTest_.is_null ()) {
581  convTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
582  }
583  }
584 
585  // Check for a change in scaling, if so we need to build new residual tests.
586  bool newResTest = false;
587  {
588  // "Residual Scaling" is the old parameter name; "Implicit
589  // Residual Scaling" is the new name. We support both options for
590  // backwards compatibility.
591  std::string tempResScale = resScale_;
592  bool implicitResidualScalingName = false;
593  if (params->isParameter ("Residual Scaling")) {
594  tempResScale = params->get<std::string> ("Residual Scaling");
595  }
596  else if (params->isParameter ("Implicit Residual Scaling")) {
597  tempResScale = params->get<std::string> ("Implicit Residual Scaling");
598  implicitResidualScalingName = true;
599  }
600 
601  // Only update the scaling if it's different.
602  if (resScale_ != tempResScale) {
603  const Belos::ScaleType resScaleType =
604  convertStringToScaleType (tempResScale);
605  resScale_ = tempResScale;
606 
607  // Update parameter in our list and residual tests, using the
608  // given parameter name.
609  if (implicitResidualScalingName) {
610  params_->set ("Implicit Residual Scaling", resScale_);
611  }
612  else {
613  params_->set ("Residual Scaling", resScale_);
614  }
615 
616  if (! convTest_.is_null ()) {
617  try {
618  convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
619  }
620  catch (std::exception& e) {
621  // Make sure the convergence test gets constructed again.
622  newResTest = true;
623  }
624  }
625  }
626  }
627 
628  // Get the deflation quorum, or number of converged systems before deflation is allowed
629  if (params->isParameter ("Deflation Quorum")) {
630  defQuorum_ = params->get ("Deflation Quorum", defQuorum_);
631  params_->set ("Deflation Quorum", defQuorum_);
632  if (! convTest_.is_null ()) {
633  convTest_->setQuorum( defQuorum_ );
634  }
635  }
636 
637  // Create status tests if we need to.
638 
639  // Basic test checks maximum iterations and native residual.
640  if (maxIterTest_.is_null ()) {
641  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
642  }
643 
644  // Implicit residual test, using the native residual to determine if convergence was achieved.
645  if (convTest_.is_null () || newResTest) {
646  convTest_ = rcp (new StatusTestResNorm_t (convtol_, defQuorum_, showMaxResNormOnly_));
647  convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
648  }
649 
650  if (sTest_.is_null () || newResTest) {
651  sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
652  }
653 
654  if (outputTest_.is_null () || newResTest) {
655  // Create the status test output class.
656  // This class manages and formats the output from the status test.
657  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
658  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
660 
661  // Set the solver string for the output test
662  const std::string solverDesc = " Pseudo Block CG ";
663  outputTest_->setSolverDesc (solverDesc);
664  }
665 
666  // Create the timer if we need to.
667  if (timerSolve_.is_null ()) {
668  const std::string solveLabel =
669  label_ + ": PseudoBlockCGSolMgr total solve time";
670 #ifdef BELOS_TEUCHOS_TIME_MONITOR
671  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
672 #endif
673  }
674 
675  // Inform the solver manager that the current parameters were set.
676  isSet_ = true;
677 }
678 
679 
680 template<class ScalarType, class MV, class OP>
681 Teuchos::RCP<const Teuchos::ParameterList>
683 {
684  using Teuchos::ParameterList;
685  using Teuchos::parameterList;
686  using Teuchos::RCP;
687 
688  if (validParams_.is_null()) {
689  // Set all the valid parameters and their default values.
690  RCP<ParameterList> pl = parameterList ();
691  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
692  "The relative residual tolerance that needs to be achieved by the\n"
693  "iterative solver in order for the linear system to be declared converged.");
694  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
695  "The maximum number of block iterations allowed for each\n"
696  "set of RHS solved.");
697  pl->set("Assert Positive Definiteness", static_cast<bool>(assertPositiveDefiniteness_default_),
698  "Whether or not to assert that the linear operator\n"
699  "and the preconditioner are indeed positive definite.");
700  pl->set("Verbosity", static_cast<int>(verbosity_default_),
701  "What type(s) of solver information should be outputted\n"
702  "to the output stream.");
703  pl->set("Output Style", static_cast<int>(outputStyle_default_),
704  "What style is used for the solver information outputted\n"
705  "to the output stream.");
706  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
707  "How often convergence information should be outputted\n"
708  "to the output stream.");
709  pl->set("Deflation Quorum", static_cast<int>(defQuorum_default_),
710  "The number of linear systems that need to converge before\n"
711  "they are deflated. This number should be <= block size.");
712  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
713  "A reference-counted pointer to the output stream where all\n"
714  "solver output is sent.");
715  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
716  "When convergence information is printed, only show the maximum\n"
717  "relative residual norm when the block size is greater than one.");
718  pl->set("Implicit Residual Scaling", resScale_default_,
719  "The type of scaling used in the residual convergence test.");
720  pl->set("Estimate Condition Number", static_cast<bool>(genCondEst_default_),
721  "Whether or not to estimate the condition number of the preconditioned system.");
722  // We leave the old name as a valid parameter for backwards
723  // compatibility (so that validateParametersAndSetDefaults()
724  // doesn't raise an exception if it encounters "Residual
725  // Scaling"). The new name was added for compatibility with other
726  // solvers, none of which use "Residual Scaling".
727  pl->set("Residual Scaling", resScale_default_,
728  "The type of scaling used in the residual convergence test. This "
729  "name is deprecated; the new name is \"Implicit Residual Scaling\".");
730  pl->set("Timer Label", static_cast<const char *>(label_default_),
731  "The string to use as a prefix for the timer labels.");
732  validParams_ = pl;
733  }
734  return validParams_;
735 }
736 
737 
738 // solve()
739 template<class ScalarType, class MV, class OP>
741 {
742  const char prefix[] = "Belos::PseudoBlockCGSolMgr::solve: ";
743 
744  // Set the current parameters if they were not set before.
745  // NOTE: This may occur if the user generated the solver manager with the default constructor and
746  // then didn't set any parameters using setParameters().
747  if (!isSet_) { setParameters( params_ ); }
748 
749  Teuchos::BLAS<int,ScalarType> blas;
750 
751  TEUCHOS_TEST_FOR_EXCEPTION
752  (! problem_->isProblemSet (), PseudoBlockCGSolMgrLinearProblemFailure,
753  prefix << "The linear problem to solve is not ready. You must call "
754  "setProblem() on the Belos::LinearProblem instance before telling the "
755  "Belos solver to solve it.");
756 
757  // Create indices for the linear systems to be solved.
758  int startPtr = 0;
759  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
760  int numCurrRHS = numRHS2Solve;
761 
762  std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
763  for (int i=0; i<numRHS2Solve; ++i) {
764  currIdx[i] = startPtr+i;
765  currIdx2[i]=i;
766  }
767 
768  // Inform the linear problem of the current linear system to solve.
769  problem_->setLSIndex( currIdx );
770 
772  // Parameter list (iteration)
773  Teuchos::ParameterList plist;
774 
775  plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
776  if(genCondEst_) plist.set("Max Size For Condest",maxIters_);
777 
778  // Reset the status test.
779  outputTest_->reset();
780 
781  // Assume convergence is achieved, then let any failed convergence set this to false.
782  bool isConverged = true;
783 
785  // Pseudo-Block CG solver
786  Teuchos::RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
787  if (numRHS2Solve == 1) {
788  block_cg_iter =
789  Teuchos::rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
790  } else {
791  block_cg_iter =
792  Teuchos::rcp (new PseudoBlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
793  }
794 
795 
796  // Enter solve() iterations
797  {
798 #ifdef BELOS_TEUCHOS_TIME_MONITOR
799  Teuchos::TimeMonitor slvtimer(*timerSolve_);
800 #endif
801 
802  bool first_time=true;
803  while ( numRHS2Solve > 0 ) {
804  if(genCondEst_ && first_time) block_cg_iter->setDoCondEst(true);
805  else block_cg_iter->setDoCondEst(false);
806 
807  // Reset the active / converged vectors from this block
808  std::vector<int> convRHSIdx;
809  std::vector<int> currRHSIdx( currIdx );
810  currRHSIdx.resize(numCurrRHS);
811 
812  // Reset the number of iterations.
813  block_cg_iter->resetNumIters();
814 
815  // Reset the number of calls that the status test output knows about.
816  outputTest_->resetNumCalls();
817 
818  // Get the current residual for this block of linear systems.
819  Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
820 
821  // Get a new state struct and initialize the solver.
823  newState.R = R_0;
824  block_cg_iter->initializeCG(newState);
825 
826  while(1) {
827 
828  // tell block_gmres_iter to iterate
829  try {
830 
831  block_cg_iter->iterate();
832 
834  //
835  // check convergence first
836  //
838  if ( convTest_->getStatus() == Passed ) {
839 
840  // Figure out which linear systems converged.
841  std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
842 
843  // If the number of converged linear systems is equal to the
844  // number of current linear systems, then we are done with this block.
845  if (convIdx.size() == currRHSIdx.size())
846  break; // break from while(1){block_cg_iter->iterate()}
847 
848  // Inform the linear problem that we are finished with this current linear system.
849  problem_->setCurrLS();
850 
851  // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
852  int have = 0;
853  std::vector<int> unconvIdx(currRHSIdx.size());
854  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
855  bool found = false;
856  for (unsigned int j=0; j<convIdx.size(); ++j) {
857  if (currRHSIdx[i] == convIdx[j]) {
858  found = true;
859  break;
860  }
861  }
862  if (!found) {
863  currIdx2[have] = currIdx2[i];
864  currRHSIdx[have++] = currRHSIdx[i];
865  }
866  }
867  currRHSIdx.resize(have);
868  currIdx2.resize(have);
869 
870  // Set the remaining indices after deflation.
871  problem_->setLSIndex( currRHSIdx );
872 
873  // Get the current residual vector.
874  std::vector<MagnitudeType> norms;
875  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
876  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
877 
878  // Set the new state and initialize the solver.
880  defstate.R = R_0;
881  block_cg_iter->initializeCG(defstate);
882  }
883 
885  //
886  // check for maximum iterations
887  //
889  else if ( maxIterTest_->getStatus() == Passed ) {
890  // we don't have convergence
891  isConverged = false;
892  break; // break from while(1){block_cg_iter->iterate()}
893  }
894 
896  //
897  // we returned from iterate(), but none of our status tests Passed.
898  // something is wrong, and it is probably our fault.
899  //
901 
902  else {
903  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
904  "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
905  }
906  }
907  catch (const std::exception &e) {
908  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
909  << block_cg_iter->getNumIters() << std::endl
910  << e.what() << std::endl;
911  throw;
912  }
913  }
914 
915  // Inform the linear problem that we are finished with this block linear system.
916  problem_->setCurrLS();
917 
918  // Update indices for the linear systems to be solved.
919  startPtr += numCurrRHS;
920  numRHS2Solve -= numCurrRHS;
921 
922  if ( numRHS2Solve > 0 ) {
923 
924  numCurrRHS = numRHS2Solve;
925  currIdx.resize( numCurrRHS );
926  currIdx2.resize( numCurrRHS );
927  for (int i=0; i<numCurrRHS; ++i)
928  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
929 
930  // Set the next indices.
931  problem_->setLSIndex( currIdx );
932  }
933  else {
934  currIdx.resize( numRHS2Solve );
935  }
936 
937  first_time=false;
938  }// while ( numRHS2Solve > 0 )
939 
940  }
941 
942  // print final summary
943  sTest_->print( printer_->stream(FinalSummary) );
944 
945  // print timing information
946 #ifdef BELOS_TEUCHOS_TIME_MONITOR
947  // Calling summarize() can be expensive, so don't call unless the
948  // user wants to print out timing details. summarize() will do all
949  // the work even if it's passed a "black hole" output stream.
950  if (verbosity_ & TimingDetails)
951  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
952 #endif
953 
954  // get iteration information for this solve
955  numIters_ = maxIterTest_->getNumIters();
956 
957  // Save the convergence test value ("achieved tolerance") for this
958  // solve.
959  const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
960  if (pTestValues != NULL && pTestValues->size () > 0) {
961  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
962  }
963 
964  // Do condition estimate, if needed
965  if (genCondEst_) {
966  ScalarType l_min, l_max;
967  Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
968  Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
969  compute_condnum_tridiag_sym(diag,offdiag,l_min,l_max,condEstimate_);
970  }
971 
972  if (! isConverged) {
973  return Unconverged; // return from PseudoBlockCGSolMgr::solve()
974  }
975  return Converged; // return from PseudoBlockCGSolMgr::solve()
976 }
977 
978 // This method requires the solver manager to return a std::string that describes itself.
979 template<class ScalarType, class MV, class OP>
981 {
982  std::ostringstream oss;
983  oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
984  oss << "{";
985  oss << "}";
986  return oss.str();
987 }
988 
989 
990 template<class ScalarType, class MV, class OP>
991 void
993 compute_condnum_tridiag_sym (Teuchos::ArrayView<MagnitudeType> diag,
994  Teuchos::ArrayView<MagnitudeType> offdiag,
995  ScalarType & lambda_min,
996  ScalarType & lambda_max,
997  ScalarType & ConditionNumber )
998 {
999  typedef Teuchos::ScalarTraits<ScalarType> STS;
1000 
1001  /* Copied from az_cg.c: compute_condnum_tridiag_sym */
1002  /* diag == ScalarType vector of size N, containing the diagonal
1003  elements of A
1004  offdiag == ScalarType vector of size N-1, containing the offdiagonal
1005  elements of A. Note that A is supposed to be symmatric
1006  */
1007  int info = 0;
1008  ScalarType scalar_dummy;
1009  MagnitudeType mag_dummy;
1010  char char_N = 'N';
1011  Teuchos::LAPACK<int,ScalarType> lapack;
1012  const int N = diag.size ();
1013 
1014  lambda_min = STS::one ();
1015  lambda_max = STS::one ();
1016  if( N > 2 ) {
1017  lapack.STEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
1018  &scalar_dummy, 1, &mag_dummy, &info);
1019  TEUCHOS_TEST_FOR_EXCEPTION
1020  (info < 0, std::logic_error, "Belos::PseudoBlockCGSolMgr::"
1021  "compute_condnum_tridiag_sym: LAPACK's _STEQR failed with info = "
1022  << info << " < 0. This suggests there might be a bug in the way Belos "
1023  "is calling LAPACK. Please report this to the Belos developers.");
1024  lambda_min = Teuchos::as<ScalarType> (diag[0]);
1025  lambda_max = Teuchos::as<ScalarType> (diag[N-1]);
1026  }
1027 
1028  // info > 0 means that LAPACK's eigensolver didn't converge. This
1029  // is unlikely but may be possible. In that case, the best we can
1030  // do is use the eigenvalues that it computes, as long as lambda_max
1031  // >= lambda_min.
1032  if (STS::real (lambda_max) < STS::real (lambda_min)) {
1033  ConditionNumber = STS::one ();
1034  }
1035  else {
1036  // It's OK for the condition number to be Inf.
1037  ConditionNumber = lambda_max / lambda_min;
1038  }
1039 
1040 } /* compute_condnum_tridiag_sym */
1041 
1042 
1043 
1044 
1045 
1046 } // end Belos namespace
1047 
1048 #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:88
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:119
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:78
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.
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:205
Pure virtual base class which describes the basic interface for a solver manager. ...
static constexpr double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:294
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:154
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:286
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