Belos Package Browser (Single Doxygen Collection)  Development
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,
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:
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.
315 
316  // Status test.
321 
322  // Current parameter list.
324 
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_;
352  ScalarType condEstimate_;
353 
354  // Timers.
355  std::string label_;
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 
752  (! problem_->isProblemSet (), PseudoBlockCGSolMgrLinearProblemFailure,
753  prefix << "The linear problem to solve is not ready. You must call "
754  "setProblem() on the Belos::LinearProblem instance before telling the "
755  "Belos solver to solve it.");
756 
757  // Create indices for the linear systems to be solved.
758  int startPtr = 0;
759  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
760  int numCurrRHS = numRHS2Solve;
761 
762  std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
763  for (int i=0; i<numRHS2Solve; ++i) {
764  currIdx[i] = startPtr+i;
765  currIdx2[i]=i;
766  }
767 
768  // Inform the linear problem of the current linear system to solve.
769  problem_->setLSIndex( currIdx );
770 
772  // Parameter list (iteration)
774 
775  plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
776  if(genCondEst_) plist.set("Max Size For Condest",maxIters_);
777 
778  // Reset the status test.
779  outputTest_->reset();
780 
781  // Assume convergence is achieved, then let any failed convergence set this to false.
782  bool isConverged = true;
783 
785  // Pseudo-Block CG solver
787  if (numRHS2Solve == 1) {
788  block_cg_iter =
789  Teuchos::rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
790  } else {
791  block_cg_iter =
792  Teuchos::rcp (new PseudoBlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_, plist));
793  }
794 
795 
796  // Enter solve() iterations
797  {
798 #ifdef BELOS_TEUCHOS_TIME_MONITOR
799  Teuchos::TimeMonitor slvtimer(*timerSolve_);
800 #endif
801 
802  bool first_time=true;
803  while ( numRHS2Solve > 0 ) {
804  if(genCondEst_ && first_time) block_cg_iter->setDoCondEst(true);
805  else block_cg_iter->setDoCondEst(false);
806 
807  // Reset the active / converged vectors from this block
808  std::vector<int> convRHSIdx;
809  std::vector<int> currRHSIdx( currIdx );
810  currRHSIdx.resize(numCurrRHS);
811 
812  // Reset the number of iterations.
813  block_cg_iter->resetNumIters();
814 
815  // Reset the number of calls that the status test output knows about.
816  outputTest_->resetNumCalls();
817 
818  // Get the current residual for this block of linear systems.
819  Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
820 
821  // Get a new state struct and initialize the solver.
823  newState.R = R_0;
824  block_cg_iter->initializeCG(newState);
825 
826  while(1) {
827 
828  // tell block_gmres_iter to iterate
829  try {
830 
831  block_cg_iter->iterate();
832 
834  //
835  // check convergence first
836  //
838  if ( convTest_->getStatus() == Passed ) {
839 
840  // Figure out which linear systems converged.
841  std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
842 
843  // If the number of converged linear systems is equal to the
844  // number of current linear systems, then we are done with this block.
845  if (convIdx.size() == currRHSIdx.size())
846  break; // break from while(1){block_cg_iter->iterate()}
847 
848  // Inform the linear problem that we are finished with this current linear system.
849  problem_->setCurrLS();
850 
851  // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
852  int have = 0;
853  std::vector<int> unconvIdx(currRHSIdx.size());
854  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
855  bool found = false;
856  for (unsigned int j=0; j<convIdx.size(); ++j) {
857  if (currRHSIdx[i] == convIdx[j]) {
858  found = true;
859  break;
860  }
861  }
862  if (!found) {
863  currIdx2[have] = currIdx2[i];
864  currRHSIdx[have++] = currRHSIdx[i];
865  }
866  }
867  currRHSIdx.resize(have);
868  currIdx2.resize(have);
869 
870  // Set the remaining indices after deflation.
871  problem_->setLSIndex( currRHSIdx );
872 
873  // Get the current residual vector.
874  std::vector<MagnitudeType> norms;
875  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
876  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
877 
878  // Set the new state and initialize the solver.
880  defstate.R = R_0;
881  block_cg_iter->initializeCG(defstate);
882  }
883 
885  //
886  // check for maximum iterations
887  //
889  else if ( maxIterTest_->getStatus() == Passed ) {
890  // we don't have convergence
891  isConverged = false;
892  break; // break from while(1){block_cg_iter->iterate()}
893  }
894 
896  //
897  // we returned from iterate(), but none of our status tests Passed.
898  // something is wrong, and it is probably our fault.
899  //
901 
902  else {
903  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
904  "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
905  }
906  }
907  catch (const std::exception &e) {
908  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
909  << block_cg_iter->getNumIters() << std::endl
910  << e.what() << std::endl;
911  throw;
912  }
913  }
914 
915  // Inform the linear problem that we are finished with this block linear system.
916  problem_->setCurrLS();
917 
918  // Update indices for the linear systems to be solved.
919  startPtr += numCurrRHS;
920  numRHS2Solve -= numCurrRHS;
921 
922  if ( numRHS2Solve > 0 ) {
923 
924  numCurrRHS = numRHS2Solve;
925  currIdx.resize( numCurrRHS );
926  currIdx2.resize( numCurrRHS );
927  for (int i=0; i<numCurrRHS; ++i)
928  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
929 
930  // Set the next indices.
931  problem_->setLSIndex( currIdx );
932  }
933  else {
934  currIdx.resize( numRHS2Solve );
935  }
936 
937  first_time=false;
938  }// while ( numRHS2Solve > 0 )
939 
940  }
941 
942  // print final summary
943  sTest_->print( printer_->stream(FinalSummary) );
944 
945  // print timing information
946 #ifdef BELOS_TEUCHOS_TIME_MONITOR
947  // Calling summarize() can be expensive, so don't call unless the
948  // user wants to print out timing details. summarize() will do all
949  // the work even if it's passed a "black hole" output stream.
950  if (verbosity_ & TimingDetails)
951  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
952 #endif
953 
954  // get iteration information for this solve
955  numIters_ = maxIterTest_->getNumIters();
956 
957  // Save the convergence test value ("achieved tolerance") for this
958  // solve.
959  const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
960  if (pTestValues != NULL && pTestValues->size () > 0) {
961  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
962  }
963 
964  // Do condition estimate, if needed
965  if (genCondEst_) {
966  ScalarType l_min, l_max;
967  Teuchos::ArrayView<MagnitudeType> diag = block_cg_iter->getDiag();
968  Teuchos::ArrayView<MagnitudeType> offdiag = block_cg_iter->getOffDiag();
969  compute_condnum_tridiag_sym(diag,offdiag,l_min,l_max,condEstimate_);
970  }
971 
972  if (! isConverged) {
973  return Unconverged; // return from PseudoBlockCGSolMgr::solve()
974  }
975  return Converged; // return from PseudoBlockCGSolMgr::solve()
976 }
977 
978 // This method requires the solver manager to return a std::string that describes itself.
979 template<class ScalarType, class MV, class OP>
981 {
982  std::ostringstream oss;
983  oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
984  oss << "{";
985  oss << "}";
986  return oss.str();
987 }
988 
989 
990 template<class ScalarType, class MV, class OP>
991 void
995  ScalarType & lambda_min,
996  ScalarType & lambda_max,
997  ScalarType & ConditionNumber )
998 {
1000 
1001  /* Copied from az_cg.c: compute_condnum_tridiag_sym */
1002  /* diag == ScalarType vector of size N, containing the diagonal
1003  elements of A
1004  offdiag == ScalarType vector of size N-1, containing the offdiagonal
1005  elements of A. Note that A is supposed to be symmatric
1006  */
1007  int info = 0;
1008  ScalarType scalar_dummy;
1009  MagnitudeType mag_dummy;
1010  char char_N = 'N';
1012  const int N = diag.size ();
1013 
1014  lambda_min = STS::one ();
1015  lambda_max = STS::one ();
1016  if( N > 2 ) {
1017  lapack.STEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
1018  &scalar_dummy, 1, &mag_dummy, &info);
1020  (info < 0, std::logic_error, "Belos::PseudoBlockCGSolMgr::"
1021  "compute_condnum_tridiag_sym: LAPACK's _STEQR failed with info = "
1022  << info << " < 0. This suggests there might be a bug in the way Belos "
1023  "is calling LAPACK. Please report this to the Belos developers.");
1024  lambda_min = Teuchos::as<ScalarType> (diag[0]);
1025  lambda_max = Teuchos::as<ScalarType> (diag[N-1]);
1026  }
1027 
1028  // info > 0 means that LAPACK's eigensolver didn't converge. This
1029  // is unlikely but may be possible. In that case, the best we can
1030  // do is use the eigenvalues that it computes, as long as lambda_max
1031  // >= lambda_min.
1032  if (STS::real (lambda_max) < STS::real (lambda_min)) {
1033  ConditionNumber = STS::one ();
1034  }
1035  else {
1036  // It's OK for the condition number to be Inf.
1037  ConditionNumber = lambda_max / lambda_min;
1038  }
1039 
1040 } /* compute_condnum_tridiag_sym */
1041 
1042 
1043 
1044 
1045 
1046 } // end Belos namespace
1047 
1048 #endif /* BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP */
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:88
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
ScalarType getConditionEstimate() const
Gets the estimated condition number.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:119
PseudoBlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:78
T & get(ParameterList &l, const std::string &name)
Structure to contain pointers to CGIteration state variables.
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)
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void STEQR(const char &COMPZ, const OrdinalType &n, ScalarType *D, ScalarType *E, ScalarType *Z, const OrdinalType &ldz, ScalarType *WORK, OrdinalType *info) const
A factory class for generating StatusTestOutput objects.
PseudoBlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
An implementation of StatusTestResNorm using a family of residual norms.
size_type size() const
Belos::StatusTest class for specifying a maximum number of iterations.
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
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 validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:205
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. ...
static constexpr double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:294
Belos concrete class for performing a single-reduction conjugate-gradient (CG) iteration.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
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)
bool is_null(const RCP< T > &p)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< const Teuchos::ParameterList > validParams_
List of valid parameters and their default values.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
PseudoBlockCGSolMgrOrthoFailure(const std::string &what_arg)
bool isType(const std::string &name) const
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:154
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
Return the residual status test.
T * getRawPtr() 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...
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Teuchos::RCP< OutputManager< ScalarType > > printer_
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
PseudoBlockCGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate or...
Class which defines basic traits for the operator type.
Belos concrete class for performing the pseudo-block CG iteration.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:286
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.