Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPseudoBlockCGSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
11 #define BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 
20 #include "BelosLinearProblem.hpp"
21 #include "BelosSolverManager.hpp"
22 
24 #include "BelosCGSingleRedIter.hpp"
25 #include "BelosCGIter.hpp"
28 #include "BelosStatusTestCombo.hpp"
30 #include "BelosOutputManager.hpp"
31 #include "Teuchos_LAPACK.hpp"
32 #ifdef BELOS_TEUCHOS_TIME_MONITOR
33 #include "Teuchos_TimeMonitor.hpp"
34 #endif
35 
55 namespace Belos {
56 
58 
59 
67  PseudoBlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
68  {}};
69 
70 
71  // Partial specialization for unsupported ScalarType types.
72  // This contains a stub implementation.
73  template<class ScalarType, class MV, class OP,
74  const bool supportsScalarType =
77  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
78  Belos::Details::LapackSupportsScalar<ScalarType>::value>
79  {
80  static const bool scalarTypeIsSupported =
82  typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
83  scalarTypeIsSupported> base_type;
84 
85  public:
87  base_type ()
88  {}
91  base_type ()
92  {}
93  virtual ~PseudoBlockCGSolMgr () {}
94 
96  getResidualStatusTest() const { return Teuchos::null; }
97  };
98 
99 
100  template<class ScalarType, class MV, class OP>
101  class PseudoBlockCGSolMgr<ScalarType, MV, OP, true> :
102  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
103  {
104  private:
108  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
110 
111  public:
112 
114 
115 
122 
140 
142  virtual ~PseudoBlockCGSolMgr() {};
143 
147  }
149 
151 
152 
153  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
154  return *problem_;
155  }
156 
159  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
160 
164 
171  return Teuchos::tuple(timerSolve_);
172  }
173 
174 
185  MagnitudeType achievedTol() const override {
186  return achievedTol_;
187  }
188 
190  int getNumIters() const override {
191  return numIters_;
192  }
193 
197  bool isLOADetected() const override { return false; }
198 
202  ScalarType getConditionEstimate() const {return condEstimate_;}
203  Teuchos::ArrayRCP<MagnitudeType> getEigenEstimates() const {return eigenEstimates_;}
204 
207  getResidualStatusTest() const { return convTest_; }
208 
210 
212 
213 
215  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
216 
218  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
219 
221 
223 
224 
228  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
230 
232 
233 
251  ReturnType solve() override;
252 
254 
257 
259  std::string description() const override;
260 
262  private:
263  // Compute the condition number estimate
264  void compute_condnum_tridiag_sym(Teuchos::ArrayView<MagnitudeType> diag,
267  ScalarType & lambda_min,
268  ScalarType & lambda_max,
269  ScalarType & ConditionNumber );
270 
271  // Linear problem.
273 
274  // Output manager.
276  Teuchos::RCP<std::ostream> outputStream_;
277 
278  // Status test.
283 
284  // Current parameter list.
286 
292  mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
293 
294  // Default solver values.
295  static constexpr int maxIters_default_ = 1000;
296  static constexpr bool assertPositiveDefiniteness_default_ = true;
297  static constexpr bool showMaxResNormOnly_default_ = false;
298  static constexpr int verbosity_default_ = Belos::Errors;
299  static constexpr int outputStyle_default_ = Belos::General;
300  static constexpr int outputFreq_default_ = -1;
301  static constexpr int defQuorum_default_ = 1;
302  static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ = false;
303  static constexpr const char * resScale_default_ = "Norm of Initial Residual";
304  static constexpr const char * label_default_ = "Belos";
305  static constexpr bool genCondEst_default_ = false;
306 
307  // Current solver values.
308  MagnitudeType convtol_,achievedTol_;
309  int maxIters_, numIters_;
310  int verbosity_, outputStyle_, outputFreq_, defQuorum_;
311  bool assertPositiveDefiniteness_, showMaxResNormOnly_;
312  bool foldConvergenceDetectionIntoAllreduce_;
313  std::string resScale_;
314  bool genCondEst_;
315  ScalarType condEstimate_;
316  Teuchos::ArrayRCP<MagnitudeType> eigenEstimates_;
317 
318  // Timers.
319  std::string label_;
320  Teuchos::RCP<Teuchos::Time> timerSolve_;
321 
322  // Internal state variables.
323  bool isSet_;
324  };
325 
326 
327 // Empty Constructor
328 template<class ScalarType, class MV, class OP>
330  outputStream_(Teuchos::rcpFromRef(std::cout)),
331  convtol_(DefaultSolverParameters::convTol),
332  maxIters_(maxIters_default_),
333  numIters_(0),
334  verbosity_(verbosity_default_),
335  outputStyle_(outputStyle_default_),
336  outputFreq_(outputFreq_default_),
337  defQuorum_(defQuorum_default_),
338  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
339  showMaxResNormOnly_(showMaxResNormOnly_default_),
340  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
341  resScale_(resScale_default_),
342  genCondEst_(genCondEst_default_),
343  condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
344  label_(label_default_),
345  isSet_(false)
346 {}
347 
348 // Basic Constructor
349 template<class ScalarType, class MV, class OP>
353  problem_(problem),
354  outputStream_(Teuchos::rcpFromRef(std::cout)),
355  convtol_(DefaultSolverParameters::convTol),
356  maxIters_(maxIters_default_),
357  numIters_(0),
358  verbosity_(verbosity_default_),
359  outputStyle_(outputStyle_default_),
360  outputFreq_(outputFreq_default_),
361  defQuorum_(defQuorum_default_),
362  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
363  showMaxResNormOnly_(showMaxResNormOnly_default_),
364  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
365  resScale_(resScale_default_),
366  genCondEst_(genCondEst_default_),
367  condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
368  label_(label_default_),
369  isSet_(false)
370 {
372  problem_.is_null (), std::invalid_argument,
373  "Belos::PseudoBlockCGSolMgr two-argument constructor: "
374  "'problem' is null. You must supply a non-null Belos::LinearProblem "
375  "instance when calling this constructor.");
376 
377  if (! pl.is_null ()) {
378  // Set the parameters using the list that was passed in.
379  setParameters (pl);
380  }
381 }
382 
383 template<class ScalarType, class MV, class OP>
384 void
387 {
389  using Teuchos::parameterList;
390  using Teuchos::RCP;
391  using Teuchos::rcp;
392 
393  RCP<const ParameterList> defaultParams = this->getValidParameters ();
394 
395  // Create the internal parameter list if one doesn't already exist.
396  // Belos' solvers treat the input ParameterList to setParameters as
397  // a "delta" -- that is, a change from the current state -- so the
398  // default parameter list (if the input is null) should be empty.
399  // This explains also why Belos' solvers copy parameters one by one
400  // from the input list to the current list.
401  //
402  // Belos obfuscates the latter, because it takes the input parameter
403  // list by RCP, rather than by (nonconst) reference. The latter
404  // would make more sense, given that it doesn't actually keep the
405  // input parameter list.
406  //
407  // Note, however, that Belos still correctly triggers the "used"
408  // field of each parameter in the input list. While isParameter()
409  // doesn't (apparently) trigger the "used" flag, get() certainly
410  // does.
411 
412  if (params_.is_null ()) {
413  // Create an empty list with the same name as the default list.
414  params_ = parameterList (defaultParams->name ());
415  } else {
416  params->validateParameters (*defaultParams);
417  }
418 
419  // Check for maximum number of iterations
420  if (params->isParameter ("Maximum Iterations")) {
421  maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
422 
423  // Update parameter in our list and in status test.
424  params_->set ("Maximum Iterations", maxIters_);
425  if (! maxIterTest_.is_null ()) {
426  maxIterTest_->setMaxIters (maxIters_);
427  }
428  }
429 
430  // Check if positive definiteness assertions are to be performed
431  if (params->isParameter ("Assert Positive Definiteness")) {
432  assertPositiveDefiniteness_ =
433  params->get ("Assert Positive Definiteness",
434  assertPositiveDefiniteness_default_);
435 
436  // Update parameter in our list.
437  params_->set ("Assert Positive Definiteness", assertPositiveDefiniteness_);
438  }
439 
440  if (params->isParameter("Fold Convergence Detection Into Allreduce")) {
441  foldConvergenceDetectionIntoAllreduce_ = params->get("Fold Convergence Detection Into Allreduce",
442  foldConvergenceDetectionIntoAllreduce_default_);
443  }
444 
445  // Check to see if the timer label changed.
446  if (params->isParameter ("Timer Label")) {
447  const std::string tempLabel = params->get ("Timer Label", label_default_);
448 
449  // Update parameter in our list and solver timer
450  if (tempLabel != label_) {
451  label_ = tempLabel;
452  params_->set ("Timer Label", label_);
453  const std::string solveLabel =
454  label_ + ": PseudoBlockCGSolMgr total solve time";
455 #ifdef BELOS_TEUCHOS_TIME_MONITOR
456  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
457 #endif
458  }
459  }
460 
461  // Check for a change in verbosity level
462  if (params->isParameter ("Verbosity")) {
463  if (Teuchos::isParameterType<int> (*params, "Verbosity")) {
464  verbosity_ = params->get ("Verbosity", verbosity_default_);
465  } else {
466  verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params, "Verbosity");
467  }
468 
469  // Update parameter in our list.
470  params_->set ("Verbosity", verbosity_);
471  if (! printer_.is_null ()) {
472  printer_->setVerbosity (verbosity_);
473  }
474  }
475 
476  // Check for a change in output style
477  if (params->isParameter ("Output Style")) {
478  if (Teuchos::isParameterType<int> (*params, "Output Style")) {
479  outputStyle_ = params->get ("Output Style", outputStyle_default_);
480  } else {
481  // FIXME (mfh 29 Jul 2015) What if the type is wrong?
482  outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params, "Output Style");
483  }
484 
485  // Reconstruct the convergence test if the explicit residual test
486  // is not being used.
487  params_->set ("Output Style", outputStyle_);
488  outputTest_ = Teuchos::null;
489  }
490 
491  // output stream
492  if (params->isParameter ("Output Stream")) {
493  outputStream_ = params->get<RCP<std::ostream> > ("Output Stream");
494 
495  // Update parameter in our list.
496  params_->set ("Output Stream", outputStream_);
497  if (! printer_.is_null ()) {
498  printer_->setOStream (outputStream_);
499  }
500  }
501 
502  // frequency level
503  if (verbosity_ & Belos::StatusTestDetails) {
504  if (params->isParameter ("Output Frequency")) {
505  outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
506  }
507 
508  // Update parameter in out list and output status test.
509  params_->set ("Output Frequency", outputFreq_);
510  if (! outputTest_.is_null ()) {
511  outputTest_->setOutputFrequency (outputFreq_);
512  }
513  }
514 
515  // Condition estimate
516  if (params->isParameter ("Estimate Condition Number")) {
517  genCondEst_ = params->get ("Estimate Condition Number", genCondEst_default_);
518  }
519 
520  // Create output manager if we need to.
521  if (printer_.is_null ()) {
522  printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
523  }
524 
525  // Convergence
526  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
527  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
528 
529  // Check for convergence tolerance
530  if (params->isParameter ("Convergence Tolerance")) {
531  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
532  convtol_ = params->get ("Convergence Tolerance",
533  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
534  }
535  else {
536  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
537  }
538 
539  // Update parameter in our list and residual tests.
540  params_->set ("Convergence Tolerance", convtol_);
541  if (! convTest_.is_null ()) {
542  convTest_->setTolerance (convtol_);
543  }
544  }
545 
546  if (params->isParameter ("Show Maximum Residual Norm Only")) {
547  showMaxResNormOnly_ = params->get<bool> ("Show Maximum Residual Norm Only");
548 
549  // Update parameter in our list and residual tests
550  params_->set ("Show Maximum Residual Norm Only", showMaxResNormOnly_);
551  if (! convTest_.is_null ()) {
552  convTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
553  }
554  }
555 
556  // Check for a change in scaling, if so we need to build new residual tests.
557  bool newResTest = false;
558  {
559  // "Residual Scaling" is the old parameter name; "Implicit
560  // Residual Scaling" is the new name. We support both options for
561  // backwards compatibility.
562  std::string tempResScale = resScale_;
563  bool implicitResidualScalingName = false;
564  if (params->isParameter ("Residual Scaling")) {
565  tempResScale = params->get<std::string> ("Residual Scaling");
566  }
567  else if (params->isParameter ("Implicit Residual Scaling")) {
568  tempResScale = params->get<std::string> ("Implicit Residual Scaling");
569  implicitResidualScalingName = true;
570  }
571 
572  // Only update the scaling if it's different.
573  if (resScale_ != tempResScale) {
574  const Belos::ScaleType resScaleType =
575  convertStringToScaleType (tempResScale);
576  resScale_ = tempResScale;
577 
578  // Update parameter in our list and residual tests, using the
579  // given parameter name.
580  if (implicitResidualScalingName) {
581  params_->set ("Implicit Residual Scaling", resScale_);
582  }
583  else {
584  params_->set ("Residual Scaling", resScale_);
585  }
586 
587  if (! convTest_.is_null ()) {
588  try {
589  convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
590  }
591  catch (std::exception& e) {
592  // Make sure the convergence test gets constructed again.
593  newResTest = true;
594  }
595  }
596  }
597  }
598 
599  // Get the deflation quorum, or number of converged systems before deflation is allowed
600  if (params->isParameter ("Deflation Quorum")) {
601  defQuorum_ = params->get ("Deflation Quorum", defQuorum_);
602  params_->set ("Deflation Quorum", defQuorum_);
603  if (! convTest_.is_null ()) {
604  convTest_->setQuorum( defQuorum_ );
605  }
606  }
607 
608  // Create status tests if we need to.
609 
610  // Basic test checks maximum iterations and native residual.
611  if (maxIterTest_.is_null ()) {
612  maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
613  }
614 
615  // Implicit residual test, using the native residual to determine if convergence was achieved.
616  if (convTest_.is_null () || newResTest) {
617  convTest_ = rcp (new StatusTestResNorm_t (convtol_, defQuorum_, showMaxResNormOnly_));
618  convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
619  }
620 
621  if (sTest_.is_null () || newResTest) {
622  sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
623  }
624 
625  if (outputTest_.is_null () || newResTest) {
626  // Create the status test output class.
627  // This class manages and formats the output from the status test.
628  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
629  outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
631 
632  // Set the solver string for the output test
633  const std::string solverDesc = " Pseudo Block CG ";
634  outputTest_->setSolverDesc (solverDesc);
635  }
636 
637  // Create the timer if we need to.
638  if (timerSolve_.is_null ()) {
639  const std::string solveLabel =
640  label_ + ": PseudoBlockCGSolMgr total solve time";
641 #ifdef BELOS_TEUCHOS_TIME_MONITOR
642  timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
643 #endif
644  }
645 
646  // Inform the solver manager that the current parameters were set.
647  isSet_ = true;
648 }
649 
650 
651 template<class ScalarType, class MV, class OP>
654 {
656  using Teuchos::parameterList;
657  using Teuchos::RCP;
658 
659  if (validParams_.is_null()) {
660  // Set all the valid parameters and their default values.
661  RCP<ParameterList> pl = parameterList ();
662  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
663  "The relative residual tolerance that needs to be achieved by the\n"
664  "iterative solver in order for the linear system to be declared converged.");
665  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
666  "The maximum number of block iterations allowed for each\n"
667  "set of RHS solved.");
668  pl->set("Assert Positive Definiteness", static_cast<bool>(assertPositiveDefiniteness_default_),
669  "Whether or not to assert that the linear operator\n"
670  "and the preconditioner are indeed positive definite.");
671  pl->set("Verbosity", static_cast<int>(verbosity_default_),
672  "What type(s) of solver information should be outputted\n"
673  "to the output stream.");
674  pl->set("Output Style", static_cast<int>(outputStyle_default_),
675  "What style is used for the solver information outputted\n"
676  "to the output stream.");
677  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
678  "How often convergence information should be outputted\n"
679  "to the output stream.");
680  pl->set("Deflation Quorum", static_cast<int>(defQuorum_default_),
681  "The number of linear systems that need to converge before\n"
682  "they are deflated. This number should be <= block size.");
683  pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
684  "A reference-counted pointer to the output stream where all\n"
685  "solver output is sent.");
686  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
687  "When convergence information is printed, only show the maximum\n"
688  "relative residual norm when the block size is greater than one.");
689  pl->set("Implicit Residual Scaling", resScale_default_,
690  "The type of scaling used in the residual convergence test.");
691  pl->set("Estimate Condition Number", static_cast<bool>(genCondEst_default_),
692  "Whether or not to estimate the condition number of the preconditioned system.");
693  // We leave the old name as a valid parameter for backwards
694  // compatibility (so that validateParametersAndSetDefaults()
695  // doesn't raise an exception if it encounters "Residual
696  // Scaling"). The new name was added for compatibility with other
697  // solvers, none of which use "Residual Scaling".
698  pl->set("Residual Scaling", resScale_default_,
699  "The type of scaling used in the residual convergence test. This "
700  "name is deprecated; the new name is \"Implicit Residual Scaling\".");
701  pl->set("Timer Label", static_cast<const char *>(label_default_),
702  "The string to use as a prefix for the timer labels.");
703  pl->set("Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
704  "Merge the allreduce for convergence detection with the one for CG.\n"
705  "This saves one all-reduce, but incurs more computation.");
706  validParams_ = pl;
707  }
708  return validParams_;
709 }
710 
711 
712 // solve()
713 template<class ScalarType, class MV, class OP>
715 {
716  const char prefix[] = "Belos::PseudoBlockCGSolMgr::solve: ";
717 
718  // Set the current parameters if they were not set before.
719  // NOTE: This may occur if the user generated the solver manager with the default constructor and
720  // then didn't set any parameters using setParameters().
721  if (!isSet_) { setParameters( params_ ); }
722 
724  (! problem_->isProblemSet (), PseudoBlockCGSolMgrLinearProblemFailure,
725  prefix << "The linear problem to solve is not ready. You must call "
726  "setProblem() on the Belos::LinearProblem instance before telling the "
727  "Belos solver to solve it.");
728 
729  // Create indices for the linear systems to be solved.
730  int startPtr = 0;
731  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
732  int numCurrRHS = numRHS2Solve;
733 
734  std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
735  for (int i=0; i<numRHS2Solve; ++i) {
736  currIdx[i] = startPtr+i;
737  currIdx2[i]=i;
738  }
739 
740  // Inform the linear problem of the current linear system to solve.
741  problem_->setLSIndex( currIdx );
742 
744  // Parameter list (iteration)
746 
747  plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
748  if(genCondEst_) plist.set("Max Size For Condest",maxIters_);
749 
750  // Reset the status test.
751  outputTest_->reset();
752 
753  // Assume convergence is achieved, then let any failed convergence set this to false.
754  bool isConverged = true;
755 
757  // Pseudo-Block CG solver
759  if (numRHS2Solve == 1) {
760  plist.set("Fold Convergence Detection Into Allreduce",
761  foldConvergenceDetectionIntoAllreduce_);
762  block_cg_iter =
763  Teuchos::rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, convTest_, plist));
764  } else {
765  block_cg_iter =
766  Teuchos::rcp (new PseudoBlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
767  }
768 
769  // Setup condition estimate
770  block_cg_iter->setDoCondEst(genCondEst_);
771  bool condEstPerf = false;
772 
773  // Enter solve() iterations
774  {
775 #ifdef BELOS_TEUCHOS_TIME_MONITOR
776  Teuchos::TimeMonitor slvtimer(*timerSolve_);
777 #endif
778 
779  while ( numRHS2Solve > 0 ) {
780 
781  // Reset the active / converged vectors from this block
782  std::vector<int> convRHSIdx;
783  std::vector<int> currRHSIdx( currIdx );
784  currRHSIdx.resize(numCurrRHS);
785 
786  // Reset the number of iterations.
787  block_cg_iter->resetNumIters();
788 
789  // Reset the number of calls that the status test output knows about.
790  outputTest_->resetNumCalls();
791 
792  // Get the current residual for this block of linear systems.
793  Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
794 
795  // Get a new state struct and initialize the solver.
797  newState.R = R_0;
798  block_cg_iter->initializeCG(newState);
799 
800  while(1) {
801 
802  // tell block_gmres_iter to iterate
803  try {
804 
805  block_cg_iter->iterate();
806 
808  //
809  // check convergence first
810  //
812  if ( convTest_->getStatus() == Passed ) {
813 
814  // Figure out which linear systems converged.
815  std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
816 
817  // If the number of converged linear systems is equal to the
818  // number of current linear systems, then we are done with this block.
819  if (convIdx.size() == currRHSIdx.size())
820  break; // break from while(1){block_cg_iter->iterate()}
821 
822  // Inform the linear problem that we are finished with this current linear system.
823  problem_->setCurrLS();
824 
825  // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
826  int have = 0;
827  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
828  bool found = false;
829  for (unsigned int j=0; j<convIdx.size(); ++j) {
830  if (currRHSIdx[i] == convIdx[j]) {
831  found = true;
832  break;
833  }
834  }
835  if (!found) {
836  currIdx2[have] = currIdx2[i];
837  currRHSIdx[have++] = currRHSIdx[i];
838  }
839  }
840  currRHSIdx.resize(have);
841  currIdx2.resize(have);
842 
843  // Compute condition estimate if the very first linear system in the block has converged.
844  if (currRHSIdx[0] != 0 && genCondEst_ && !condEstPerf)
845  {
846  // Compute the estimate.
847  ScalarType l_min, l_max;
848  Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
849  Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
850  compute_condnum_tridiag_sym(diag,offdiag,eigenEstimates_,l_min,l_max,condEstimate_);
851 
852  // Make sure not to do more condition estimate computations for this solve.
853  block_cg_iter->setDoCondEst(false);
854  condEstPerf = true;
855  }
856 
857  // Set the remaining indices after deflation.
858  problem_->setLSIndex( currRHSIdx );
859 
860  // Get the current residual vector.
861  std::vector<MagnitudeType> norms;
862  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
863  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
864 
865  // Set the new state and initialize the solver.
867  defstate.R = R_0;
868  block_cg_iter->initializeCG(defstate);
869  }
870 
872  //
873  // check for maximum iterations
874  //
876  else if ( maxIterTest_->getStatus() == Passed ) {
877  // we don't have convergence
878  isConverged = false;
879  break; // break from while(1){block_cg_iter->iterate()}
880  }
881 
883  //
884  // we returned from iterate(), but none of our status tests Passed.
885  // something is wrong, and it is probably our fault.
886  //
888 
889  else {
890  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
891  "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
892  }
893  }
894  catch (const StatusTestNaNError& e) {
895  // A NaN was detected in the solver. Set the solution to zero and return unconverged.
896  achievedTol_ = MT::one();
897  Teuchos::RCP<MV> X = problem_->getLHS();
898  MVT::MvInit( *X, SCT::zero() );
899  printer_->stream(Warnings) << "Belos::PseudoBlockCGSolMgr::solve(): Warning! NaN has been detected!"
900  << std::endl;
901  return Unconverged;
902  }
903  catch (const std::exception &e) {
904  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
905  << block_cg_iter->getNumIters() << std::endl
906  << e.what() << std::endl;
907  throw;
908  }
909  }
910 
911  // Inform the linear problem that we are finished with this block linear system.
912  problem_->setCurrLS();
913 
914  // Update indices for the linear systems to be solved.
915  startPtr += numCurrRHS;
916  numRHS2Solve -= numCurrRHS;
917 
918  if ( numRHS2Solve > 0 ) {
919 
920  numCurrRHS = numRHS2Solve;
921  currIdx.resize( numCurrRHS );
922  currIdx2.resize( numCurrRHS );
923  for (int i=0; i<numCurrRHS; ++i)
924  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
925 
926  // Set the next indices.
927  problem_->setLSIndex( currIdx );
928  }
929  else {
930  currIdx.resize( numRHS2Solve );
931  }
932 
933  }// while ( numRHS2Solve > 0 )
934 
935  }
936 
937  // print final summary
938  sTest_->print( printer_->stream(FinalSummary) );
939 
940  // print timing information
941 #ifdef BELOS_TEUCHOS_TIME_MONITOR
942  // Calling summarize() can be expensive, so don't call unless the
943  // user wants to print out timing details. summarize() will do all
944  // the work even if it's passed a "black hole" output stream.
945  if (verbosity_ & TimingDetails)
946  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
947 #endif
948 
949  // get iteration information for this solve
950  numIters_ = maxIterTest_->getNumIters();
951 
952  // Save the convergence test value ("achieved tolerance") for this
953  // solve.
954  const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
955  if (pTestValues != NULL && pTestValues->size () > 0) {
956  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
957  }
958 
959  // Do condition estimate, if needed
960  if (genCondEst_ && !condEstPerf) {
961  ScalarType l_min, l_max;
962  Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
963  Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
964  compute_condnum_tridiag_sym(diag,offdiag,eigenEstimates_,l_min,l_max,condEstimate_);
965  condEstPerf = true;
966  }
967 
968  if (! isConverged) {
969  return Unconverged; // return from PseudoBlockCGSolMgr::solve()
970  }
971  return Converged; // return from PseudoBlockCGSolMgr::solve()
972 }
973 
974 // This method requires the solver manager to return a std::string that describes itself.
975 template<class ScalarType, class MV, class OP>
977 {
978  std::ostringstream oss;
979  oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
980  oss << "{";
981  oss << "}";
982  return oss.str();
983 }
984 
985 
986 template<class ScalarType, class MV, class OP>
987 void
992  ScalarType & lambda_min,
993  ScalarType & lambda_max,
994  ScalarType & ConditionNumber )
995 {
997 
998  /* Copied from az_cg.c: compute_condnum_tridiag_sym */
999  /* diag == ScalarType vector of size N, containing the diagonal
1000  elements of A
1001  offdiag == ScalarType vector of size N-1, containing the offdiagonal
1002  elements of A. Note that A is supposed to be symmatric
1003  */
1004  int info = 0;
1005  const int N = diag.size ();
1006  ScalarType scalar_dummy;
1007  std::vector<MagnitudeType> mag_dummy(4*N);
1008  char char_N = 'N';
1010 
1011  lambdas.resize(N, 0.0);
1012  lambda_min = STS::one ();
1013  lambda_max = STS::one ();
1014  if( N > 2 ) {
1015  lapack.PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
1016  &scalar_dummy, 1, &mag_dummy[0], &info);
1018  (info < 0, std::logic_error, "Belos::PseudoBlockCGSolMgr::"
1019  "compute_condnum_tridiag_sym: LAPACK's _PTEQR failed with info = "
1020  << info << " < 0. This suggests there might be a bug in the way Belos "
1021  "is calling LAPACK. Please report this to the Belos developers.");
1022  for (int k = 0; k < N; k++) {
1023  lambdas[k] = diag[N - 1 - k];
1024  }
1025  lambda_min = Teuchos::as<ScalarType> (diag[N-1]);
1026  lambda_max = Teuchos::as<ScalarType> (diag[0]);
1027  }
1028 
1029  // info > 0 means that LAPACK's eigensolver didn't converge. This
1030  // is unlikely but may be possible. In that case, the best we can
1031  // do is use the eigenvalues that it computes, as long as lambda_max
1032  // >= lambda_min.
1033  if (STS::real (lambda_max) < STS::real (lambda_min)) {
1034  ConditionNumber = STS::one ();
1035  }
1036  else {
1037  // It's OK for the condition number to be Inf.
1038  ConditionNumber = lambda_max / lambda_min;
1039  }
1040 
1041 } /* compute_condnum_tridiag_sym */
1042 
1043 
1044 
1045 
1046 
1047 } // end Belos namespace
1048 
1049 #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:74
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
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:88
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:47
Structure to contain pointers to CGIteration state variables.
T & get(const std::string &name, T def_value)
Belos concrete class for performing the conjugate-gradient (CG) iteration.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
PseudoBlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
size_type size() const
An implementation of StatusTestResNorm using a family of residual norms.
ScalarType getConditionEstimate() const
Gets the estimated condition number.
void PTEQR(const char &COMPZ, const OrdinalType &n, MagnitudeType *D, MagnitudeType *E, ScalarType *Z, const OrdinalType &ldz, MagnitudeType *WORK, OrdinalType *info) const
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:261
Belos::StatusTest class for specifying a maximum number of iterations.
static std::string name()
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.
void resize(const size_type n, const T &val=T())
bool isParameter(const std::string &name) const
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:174
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
Teuchos::ArrayRCP< MagnitudeType > getEigenEstimates() const
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
Return the residual status test.
A linear system to solve, and its associated information.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Class which describes the linear problem to be solved by the iterative solver.
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:123
T * getRawPtr() const
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
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...
bool isType(const std::string &name) const
A class for extending the status testing capabilities of Belos via logical combinations.
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:28
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:251
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.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
bool is_null() const

Generated on Tue Jul 23 2024 09:24:22 for Belos by doxygen 1.8.5