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 //
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  {}
131  base_type ()
132  {}
133  virtual ~PseudoBlockCGSolMgr () {}
134 
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:
148  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
150 
151  public:
152 
154 
155 
162 
180 
182  virtual ~PseudoBlockCGSolMgr() {};
183 
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 
204 
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 
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,
305  ScalarType & lambda_min,
306  ScalarType & lambda_max,
307  ScalarType & ConditionNumber );
308 
309  // Linear problem.
311 
312  // Output manager.
314  Teuchos::RCP<std::ostream> outputStream_;
315 
316  // Status test.
321 
322  // Current parameter list.
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>
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 {
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
421 {
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>
683 {
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 
750  (! problem_->isProblemSet (), PseudoBlockCGSolMgrLinearProblemFailure,
751  prefix << "The linear problem to solve is not ready. You must call "
752  "setProblem() on the Belos::LinearProblem instance before telling the "
753  "Belos solver to solve it.");
754 
755  // Create indices for the linear systems to be solved.
756  int startPtr = 0;
757  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
758  int numCurrRHS = numRHS2Solve;
759 
760  std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
761  for (int i=0; i<numRHS2Solve; ++i) {
762  currIdx[i] = startPtr+i;
763  currIdx2[i]=i;
764  }
765 
766  // Inform the linear problem of the current linear system to solve.
767  problem_->setLSIndex( currIdx );
768 
770  // Parameter list (iteration)
772 
773  plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
774  if(genCondEst_) plist.set("Max Size For Condest",maxIters_);
775 
776  // Reset the status test.
777  outputTest_->reset();
778 
779  // Assume convergence is achieved, then let any failed convergence set this to false.
780  bool isConverged = true;
781 
783  // Pseudo-Block CG solver
785  if (numRHS2Solve == 1) {
786  block_cg_iter =
787  Teuchos::rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
788  } else {
789  block_cg_iter =
790  Teuchos::rcp (new PseudoBlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
791  }
792 
793  // Setup condition estimate
794  block_cg_iter->setDoCondEst(genCondEst_);
795  bool condEstPerf = false;
796 
797  // Enter solve() iterations
798  {
799 #ifdef BELOS_TEUCHOS_TIME_MONITOR
800  Teuchos::TimeMonitor slvtimer(*timerSolve_);
801 #endif
802 
803  while ( numRHS2Solve > 0 ) {
804 
805  // Reset the active / converged vectors from this block
806  std::vector<int> convRHSIdx;
807  std::vector<int> currRHSIdx( currIdx );
808  currRHSIdx.resize(numCurrRHS);
809 
810  // Reset the number of iterations.
811  block_cg_iter->resetNumIters();
812 
813  // Reset the number of calls that the status test output knows about.
814  outputTest_->resetNumCalls();
815 
816  // Get the current residual for this block of linear systems.
817  Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
818 
819  // Get a new state struct and initialize the solver.
821  newState.R = R_0;
822  block_cg_iter->initializeCG(newState);
823 
824  while(1) {
825 
826  // tell block_gmres_iter to iterate
827  try {
828 
829  block_cg_iter->iterate();
830 
832  //
833  // check convergence first
834  //
836  if ( convTest_->getStatus() == Passed ) {
837 
838  // Figure out which linear systems converged.
839  std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
840 
841  // If the number of converged linear systems is equal to the
842  // number of current linear systems, then we are done with this block.
843  if (convIdx.size() == currRHSIdx.size())
844  break; // break from while(1){block_cg_iter->iterate()}
845 
846  // Inform the linear problem that we are finished with this current linear system.
847  problem_->setCurrLS();
848 
849  // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
850  int have = 0;
851  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
852  bool found = false;
853  for (unsigned int j=0; j<convIdx.size(); ++j) {
854  if (currRHSIdx[i] == convIdx[j]) {
855  found = true;
856  break;
857  }
858  }
859  if (!found) {
860  currIdx2[have] = currIdx2[i];
861  currRHSIdx[have++] = currRHSIdx[i];
862  }
863  }
864  currRHSIdx.resize(have);
865  currIdx2.resize(have);
866 
867  // Compute condition estimate if the very first linear system in the block has converged.
868  if (currRHSIdx[0] != 0 && genCondEst_ && !condEstPerf)
869  {
870  // Compute the estimate.
871  ScalarType l_min, l_max;
872  Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
873  Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
874  compute_condnum_tridiag_sym(diag,offdiag,l_min,l_max,condEstimate_);
875 
876  // Make sure not to do more condition estimate computations for this solve.
877  block_cg_iter->setDoCondEst(false);
878  condEstPerf = true;
879  }
880 
881  // Set the remaining indices after deflation.
882  problem_->setLSIndex( currRHSIdx );
883 
884  // Get the current residual vector.
885  std::vector<MagnitudeType> norms;
886  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
887  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
888 
889  // Set the new state and initialize the solver.
891  defstate.R = R_0;
892  block_cg_iter->initializeCG(defstate);
893  }
894 
896  //
897  // check for maximum iterations
898  //
900  else if ( maxIterTest_->getStatus() == Passed ) {
901  // we don't have convergence
902  isConverged = false;
903  break; // break from while(1){block_cg_iter->iterate()}
904  }
905 
907  //
908  // we returned from iterate(), but none of our status tests Passed.
909  // something is wrong, and it is probably our fault.
910  //
912 
913  else {
914  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
915  "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
916  }
917  }
918  catch (const std::exception &e) {
919  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
920  << block_cg_iter->getNumIters() << std::endl
921  << e.what() << std::endl;
922  throw;
923  }
924  }
925 
926  // Inform the linear problem that we are finished with this block linear system.
927  problem_->setCurrLS();
928 
929  // Update indices for the linear systems to be solved.
930  startPtr += numCurrRHS;
931  numRHS2Solve -= numCurrRHS;
932 
933  if ( numRHS2Solve > 0 ) {
934 
935  numCurrRHS = numRHS2Solve;
936  currIdx.resize( numCurrRHS );
937  currIdx2.resize( numCurrRHS );
938  for (int i=0; i<numCurrRHS; ++i)
939  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
940 
941  // Set the next indices.
942  problem_->setLSIndex( currIdx );
943  }
944  else {
945  currIdx.resize( numRHS2Solve );
946  }
947 
948  }// while ( numRHS2Solve > 0 )
949 
950  }
951 
952  // print final summary
953  sTest_->print( printer_->stream(FinalSummary) );
954 
955  // print timing information
956 #ifdef BELOS_TEUCHOS_TIME_MONITOR
957  // Calling summarize() can be expensive, so don't call unless the
958  // user wants to print out timing details. summarize() will do all
959  // the work even if it's passed a "black hole" output stream.
960  if (verbosity_ & TimingDetails)
961  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
962 #endif
963 
964  // get iteration information for this solve
965  numIters_ = maxIterTest_->getNumIters();
966 
967  // Save the convergence test value ("achieved tolerance") for this
968  // solve.
969  const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
970  if (pTestValues != NULL && pTestValues->size () > 0) {
971  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
972  }
973 
974  // Do condition estimate, if needed
975  if (genCondEst_ && !condEstPerf) {
976  ScalarType l_min, l_max;
977  Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
978  Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
979  compute_condnum_tridiag_sym(diag,offdiag,l_min,l_max,condEstimate_);
980  condEstPerf = true;
981  }
982 
983  if (! isConverged) {
984  return Unconverged; // return from PseudoBlockCGSolMgr::solve()
985  }
986  return Converged; // return from PseudoBlockCGSolMgr::solve()
987 }
988 
989 // This method requires the solver manager to return a std::string that describes itself.
990 template<class ScalarType, class MV, class OP>
992 {
993  std::ostringstream oss;
994  oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
995  oss << "{";
996  oss << "}";
997  return oss.str();
998 }
999 
1000 
1001 template<class ScalarType, class MV, class OP>
1002 void
1006  ScalarType & lambda_min,
1007  ScalarType & lambda_max,
1008  ScalarType & ConditionNumber )
1009 {
1011 
1012  /* Copied from az_cg.c: compute_condnum_tridiag_sym */
1013  /* diag == ScalarType vector of size N, containing the diagonal
1014  elements of A
1015  offdiag == ScalarType vector of size N-1, containing the offdiagonal
1016  elements of A. Note that A is supposed to be symmatric
1017  */
1018  int info = 0;
1019  const int N = diag.size ();
1020  ScalarType scalar_dummy;
1021  std::vector<MagnitudeType> mag_dummy(4*N);
1022  char char_N = 'N';
1024 
1025  lambda_min = STS::one ();
1026  lambda_max = STS::one ();
1027  if( N > 2 ) {
1028  lapack.PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
1029  &scalar_dummy, 1, &mag_dummy[0], &info);
1031  (info < 0, std::logic_error, "Belos::PseudoBlockCGSolMgr::"
1032  "compute_condnum_tridiag_sym: LAPACK's _PTEQR failed with info = "
1033  << info << " < 0. This suggests there might be a bug in the way Belos "
1034  "is calling LAPACK. Please report this to the Belos developers.");
1035  lambda_min = Teuchos::as<ScalarType> (diag[N-1]);
1036  lambda_max = Teuchos::as<ScalarType> (diag[0]);
1037  }
1038 
1039  // info > 0 means that LAPACK's eigensolver didn't converge. This
1040  // is unlikely but may be possible. In that case, the best we can
1041  // do is use the eigenvalues that it computes, as long as lambda_max
1042  // >= lambda_min.
1043  if (STS::real (lambda_max) < STS::real (lambda_min)) {
1044  ConditionNumber = STS::one ();
1045  }
1046  else {
1047  // It's OK for the condition number to be Inf.
1048  ConditionNumber = lambda_max / lambda_min;
1049  }
1050 
1051 } /* compute_condnum_tridiag_sym */
1052 
1053 
1054 
1055 
1056 
1057 } // end Belos namespace
1058 
1059 #endif /* BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP */
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:120
PseudoBlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:78
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:301
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.
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:206
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::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.
PseudoBlockCGSolMgrOrthoFailure(const std::string &what_arg)
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
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.
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:291
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 Fri Aug 14 2020 10:48:34 for Belos by doxygen 1.8.5