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

Generated for Belos by doxygen 1.8.14