Belos Package Browser (Single Doxygen Collection)  Development
BelosPseudoBlockStochasticCGSolMgr.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_STOCHASTIC_CG_SOLMGR_HPP
43 #define BELOS_PSEUDO_BLOCK_STOCHASTIC_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 #ifdef BELOS_TEUCHOS_TIME_MONITOR
62 #include "Teuchos_TimeMonitor.hpp"
63 #endif
64 
74 namespace Belos {
75 
77 
78 
86  PseudoBlockStochasticCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
87  {}};
88 
89  template<class ScalarType, class MV, class OP>
90  class PseudoBlockStochasticCGSolMgr : public SolverManager<ScalarType,MV,OP> {
91 
92  private:
95  typedef Teuchos::ScalarTraits<ScalarType> SCT;
96  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
97  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
98 
99  public:
100 
102 
103 
110 
120  PseudoBlockStochasticCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
121  const Teuchos::RCP<Teuchos::ParameterList> &pl );
122 
125 
127  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
128  return Teuchos::rcp(new PseudoBlockStochasticCGSolMgr<ScalarType,MV,OP>);
129  }
131 
133 
134 
135  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
136  return *problem_;
137  }
138 
141  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
142 
145  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
146 
152  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
153  return Teuchos::tuple(timerSolve_);
154  }
155 
157  int getNumIters() const override {
158  return numIters_;
159  }
160 
164  bool isLOADetected() const override { return false; }
165 
167 
169 
170 
172  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
173 
175  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
176 
178 
180 
181 
185  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
187 
189 
190 
208  ReturnType solve() override;
209 
211 
213  Teuchos::RCP<MV> getStochasticVector() { return Y_;}
214 
217 
219  std::string description() const override;
220 
222 
223  private:
224 
225  // Linear problem.
226  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
227 
228  // Output manager.
229  Teuchos::RCP<OutputManager<ScalarType> > printer_;
230  Teuchos::RCP<std::ostream> outputStream_;
231 
232  // Status test.
233  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
234  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
235  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
236  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
237 
238  // Current parameter list.
239  Teuchos::RCP<Teuchos::ParameterList> params_;
240 
246  mutable Teuchos::RCP<const Teuchos::ParameterList> validParams_;
247 
248  // Default solver values.
249  static constexpr int maxIters_default_ = 1000;
250  static constexpr bool assertPositiveDefiniteness_default_ = true;
251  static constexpr bool showMaxResNormOnly_default_ = false;
252  static constexpr int verbosity_default_ = Belos::Errors;
253  static constexpr int outputStyle_default_ = Belos::General;
254  static constexpr int outputFreq_default_ = -1;
255  static constexpr int defQuorum_default_ = 1;
256  static constexpr const char * resScale_default_ = "Norm of Initial Residual";
257  static constexpr const char * label_default_ = "Belos";
258 // https://stackoverflow.com/questions/24398102/constexpr-and-initialization-of-a-static-const-void-pointer-with-reinterpret-cas
259 #if defined(_WIN32) && defined(__clang__)
260  static constexpr std::ostream * outputStream_default_ =
261  __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
262 #else
263  static constexpr std::ostream * outputStream_default_ = &std::cout;
264 #endif
265 
266  // Current solver values.
271  std::string resScale_;
272 
273  // Timers.
274  std::string label_;
275  Teuchos::RCP<Teuchos::Time> timerSolve_;
276 
277  // Internal state variables.
278  bool isSet_;
279 
280  // Stashed copy of the stochastic vector
281  Teuchos::RCP<MV> Y_;
282 
283  };
284 
285 
286 // Empty Constructor
287 template<class ScalarType, class MV, class OP>
289  outputStream_(Teuchos::rcp(outputStream_default_,false)),
290  convtol_(DefaultSolverParameters::convTol),
291  maxIters_(maxIters_default_),
292  numIters_(0),
293  verbosity_(verbosity_default_),
294  outputStyle_(outputStyle_default_),
295  outputFreq_(outputFreq_default_),
296  defQuorum_(defQuorum_default_),
297  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
298  showMaxResNormOnly_(showMaxResNormOnly_default_),
299  resScale_(resScale_default_),
300  label_(label_default_),
301  isSet_(false)
302 {}
303 
304 // Basic Constructor
305 template<class ScalarType, class MV, class OP>
308  const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
309  problem_(problem),
310  outputStream_(Teuchos::rcp(outputStream_default_,false)),
311  convtol_(DefaultSolverParameters::convTol),
312  maxIters_(maxIters_default_),
313  numIters_(0),
314  verbosity_(verbosity_default_),
315  outputStyle_(outputStyle_default_),
316  outputFreq_(outputFreq_default_),
317  defQuorum_(defQuorum_default_),
318  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
319  showMaxResNormOnly_(showMaxResNormOnly_default_),
320  resScale_(resScale_default_),
321  label_(label_default_),
322  isSet_(false)
323 {
324  TEUCHOS_TEST_FOR_EXCEPTION(
325  problem_.is_null (), std::invalid_argument,
326  "Belos::PseudoBlockStochasticCGSolMgr two-argument constructor: "
327  "'problem' is null. You must supply a non-null Belos::LinearProblem "
328  "instance when calling this constructor.");
329 
330  if (! pl.is_null ()) {
331  // Set the parameters using the list that was passed in.
332  setParameters (pl);
333  }
334 }
335 
336 template<class ScalarType, class MV, class OP>
337 void PseudoBlockStochasticCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
338 {
339  using Teuchos::ParameterList;
340  using Teuchos::parameterList;
341  using Teuchos::RCP;
342 
343  RCP<const ParameterList> defaultParams = getValidParameters();
344 
345  // Create the internal parameter list if one doesn't already exist.
346  if (params_.is_null()) {
347  params_ = parameterList (*defaultParams);
348  } else {
349  params->validateParameters (*defaultParams);
350  }
351 
352  // Check for maximum number of iterations
353  if (params->isParameter("Maximum Iterations")) {
354  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
355 
356  // Update parameter in our list and in status test.
357  params_->set("Maximum Iterations", maxIters_);
358  if (maxIterTest_!=Teuchos::null)
359  maxIterTest_->setMaxIters( maxIters_ );
360  }
361 
362  // Check if positive definiteness assertions are to be performed
363  if (params->isParameter("Assert Positive Definiteness")) {
364  assertPositiveDefiniteness_ = params->get("Assert Positive Definiteness",assertPositiveDefiniteness_default_);
365 
366  // Update parameter in our list.
367  params_->set("Assert Positive Definiteness", assertPositiveDefiniteness_);
368  }
369 
370  // Check to see if the timer label changed.
371  if (params->isParameter("Timer Label")) {
372  std::string tempLabel = params->get("Timer Label", label_default_);
373 
374  // Update parameter in our list and solver timer
375  if (tempLabel != label_) {
376  label_ = tempLabel;
377  params_->set("Timer Label", label_);
378  std::string solveLabel = label_ + ": PseudoBlockStochasticCGSolMgr total solve time";
379 #ifdef BELOS_TEUCHOS_TIME_MONITOR
380  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
381 #endif
382  }
383  }
384 
385  // Check for a change in verbosity level
386  if (params->isParameter("Verbosity")) {
387  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
388  verbosity_ = params->get("Verbosity", verbosity_default_);
389  } else {
390  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
391  }
392 
393  // Update parameter in our list.
394  params_->set("Verbosity", verbosity_);
395  if (printer_ != Teuchos::null)
396  printer_->setVerbosity(verbosity_);
397  }
398 
399  // Check for a change in output style
400  if (params->isParameter("Output Style")) {
401  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
402  outputStyle_ = params->get("Output Style", outputStyle_default_);
403  } else {
404  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
405  }
406 
407  // Reconstruct the convergence test if the explicit residual test is not being used.
408  params_->set("Output Style", outputStyle_);
409  outputTest_ = Teuchos::null;
410  }
411 
412  // output stream
413  if (params->isParameter("Output Stream")) {
414  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
415 
416  // Update parameter in our list.
417  params_->set("Output Stream", outputStream_);
418  if (printer_ != Teuchos::null)
419  printer_->setOStream( outputStream_ );
420  }
421 
422  // frequency level
423  if (verbosity_ & Belos::StatusTestDetails) {
424  if (params->isParameter("Output Frequency")) {
425  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
426  }
427 
428  // Update parameter in out list and output status test.
429  params_->set("Output Frequency", outputFreq_);
430  if (outputTest_ != Teuchos::null)
431  outputTest_->setOutputFrequency( outputFreq_ );
432  }
433 
434  // Create output manager if we need to.
435  if (printer_ == Teuchos::null) {
436  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
437  }
438 
439  // Convergence
440  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
441  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
442 
443  // Check for convergence tolerance
444  if (params->isParameter("Convergence Tolerance")) {
445  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
446  convtol_ = params->get ("Convergence Tolerance",
447  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
448  }
449  else {
450  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
451  }
452 
453  // Update parameter in our list and residual tests.
454  params_->set("Convergence Tolerance", convtol_);
455  if (convTest_ != Teuchos::null)
456  convTest_->setTolerance( convtol_ );
457  }
458 
459  if (params->isParameter("Show Maximum Residual Norm Only")) {
460  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
461 
462  // Update parameter in our list and residual tests
463  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
464  if (convTest_ != Teuchos::null)
465  convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
466  }
467 
468  // Check for a change in scaling, if so we need to build new residual tests.
469  bool newResTest = false;
470  {
471  // "Residual Scaling" is the old parameter name; "Implicit
472  // Residual Scaling" is the new name. We support both options for
473  // backwards compatibility.
474  std::string tempResScale = resScale_;
475  bool implicitResidualScalingName = false;
476  if (params->isParameter ("Residual Scaling")) {
477  tempResScale = params->get<std::string> ("Residual Scaling");
478  }
479  else if (params->isParameter ("Implicit Residual Scaling")) {
480  tempResScale = params->get<std::string> ("Implicit Residual Scaling");
481  implicitResidualScalingName = true;
482  }
483 
484  // Only update the scaling if it's different.
485  if (resScale_ != tempResScale) {
486  Belos::ScaleType resScaleType = convertStringToScaleType( tempResScale );
487  resScale_ = tempResScale;
488 
489  // Update parameter in our list and residual tests, using the
490  // given parameter name.
491  if (implicitResidualScalingName) {
492  params_->set ("Implicit Residual Scaling", resScale_);
493  }
494  else {
495  params_->set ("Residual Scaling", resScale_);
496  }
497 
498  if (! convTest_.is_null()) {
499  try {
500  convTest_->defineScaleForm( resScaleType, Belos::TwoNorm );
501  }
502  catch (std::exception& e) {
503  // Make sure the convergence test gets constructed again.
504  newResTest = true;
505  }
506  }
507  }
508  }
509 
510  // Get the deflation quorum, or number of converged systems before deflation is allowed
511  if (params->isParameter("Deflation Quorum")) {
512  defQuorum_ = params->get("Deflation Quorum", defQuorum_);
513  params_->set("Deflation Quorum", defQuorum_);
514  if (convTest_ != Teuchos::null)
515  convTest_->setQuorum( defQuorum_ );
516  }
517 
518  // Create status tests if we need to.
519 
520  // Basic test checks maximum iterations and native residual.
521  if (maxIterTest_ == Teuchos::null)
522  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
523 
524  // Implicit residual test, using the native residual to determine if convergence was achieved.
525  if (convTest_ == Teuchos::null || newResTest) {
526  convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, defQuorum_, showMaxResNormOnly_ ) );
527  convTest_->defineScaleForm( convertStringToScaleType( resScale_ ), Belos::TwoNorm );
528  }
529 
530  if (sTest_ == Teuchos::null || newResTest)
531  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
532 
533  if (outputTest_ == Teuchos::null || newResTest) {
534 
535  // Create the status test output class.
536  // This class manages and formats the output from the status test.
537  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
538  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
539 
540  // Set the solver string for the output test
541  std::string solverDesc = " Pseudo Block CG ";
542  outputTest_->setSolverDesc( solverDesc );
543 
544  }
545 
546  // Create the timer if we need to.
547  if (timerSolve_ == Teuchos::null) {
548  std::string solveLabel = label_ + ": PseudoBlockStochasticCGSolMgr total solve time";
549 #ifdef BELOS_TEUCHOS_TIME_MONITOR
550  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
551 #endif
552  }
553 
554  // Inform the solver manager that the current parameters were set.
555  isSet_ = true;
556 }
557 
558 
559 template<class ScalarType, class MV, class OP>
560 Teuchos::RCP<const Teuchos::ParameterList>
562 {
563  using Teuchos::ParameterList;
564  using Teuchos::parameterList;
565  using Teuchos::RCP;
566 
567  if (validParams_.is_null()) {
568  // Set all the valid parameters and their default values.
569 
570  // The static_cast is to resolve an issue with older clang versions which
571  // would cause the constexpr to link fail. With c++17 the problem is resolved.
572  RCP<ParameterList> pl = parameterList ();
573  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
574  "The relative residual tolerance that needs to be achieved by the\n"
575  "iterative solver in order for the linera system to be declared converged.");
576  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
577  "The maximum number of block iterations allowed for each\n"
578  "set of RHS solved.");
579  pl->set("Assert Positive Definiteness", static_cast<bool>(assertPositiveDefiniteness_default_),
580  "Whether or not to assert that the linear operator\n"
581  "and the preconditioner are indeed positive definite.");
582  pl->set("Verbosity", static_cast<int>(verbosity_default_),
583  "What type(s) of solver information should be outputted\n"
584  "to the output stream.");
585  pl->set("Output Style", static_cast<int>(outputStyle_default_),
586  "What style is used for the solver information outputted\n"
587  "to the output stream.");
588  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
589  "How often convergence information should be outputted\n"
590  "to the output stream.");
591  pl->set("Deflation Quorum", static_cast<int>(defQuorum_default_),
592  "The number of linear systems that need to converge before\n"
593  "they are deflated. This number should be <= block size.");
594  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
595  "A reference-counted pointer to the output stream where all\n"
596  "solver output is sent.");
597  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
598  "When convergence information is printed, only show the maximum\n"
599  "relative residual norm when the block size is greater than one.");
600  pl->set("Implicit Residual Scaling", resScale_default_,
601  "The type of scaling used in the residual convergence test.");
602  // We leave the old name as a valid parameter for backwards
603  // compatibility (so that validateParametersAndSetDefaults()
604  // doesn't raise an exception if it encounters "Residual
605  // Scaling"). The new name was added for compatibility with other
606  // solvers, none of which use "Residual Scaling".
607  pl->set("Residual Scaling", resScale_default_,
608  "The type of scaling used in the residual convergence test. This "
609  "name is deprecated; the new name is \"Implicit Residual Scaling\".");
610  pl->set("Timer Label", static_cast<const char *>(label_default_),
611  "The string to use as a prefix for the timer labels.");
612  validParams_ = pl;
613  }
614  return validParams_;
615 }
616 
617 
618 // solve()
619 template<class ScalarType, class MV, class OP>
621 
622  // Set the current parameters if they were not set before.
623  // NOTE: This may occur if the user generated the solver manager with the default constructor and
624  // then didn't set any parameters using setParameters().
625  if (!isSet_) { setParameters( params_ ); }
626 
627  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockStochasticCGSolMgrLinearProblemFailure,
628  "Belos::PseudoBlockStochasticCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
629 
630  // Create indices for the linear systems to be solved.
631  int startPtr = 0;
632  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
633  int numCurrRHS = numRHS2Solve;
634 
635  std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
636  for (int i=0; i<numRHS2Solve; ++i) {
637  currIdx[i] = startPtr+i;
638  currIdx2[i]=i;
639  }
640 
641  // Inform the linear problem of the current linear system to solve.
642  problem_->setLSIndex( currIdx );
643 
645  // Parameter list
646  Teuchos::ParameterList plist;
647 
648  plist.set("Assert Positive Definiteness",assertPositiveDefiniteness_);
649 
650  // Reset the status test.
651  outputTest_->reset();
652 
653  // Assume convergence is achieved, then let any failed convergence set this to false.
654  bool isConverged = true;
655 
657  // Pseudo-Block CG solver
658 
659  Teuchos::RCP<PseudoBlockStochasticCGIter<ScalarType,MV,OP> > block_cg_iter
660  = Teuchos::rcp( new PseudoBlockStochasticCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
661 
662  // Enter solve() iterations
663  {
664 #ifdef BELOS_TEUCHOS_TIME_MONITOR
665  Teuchos::TimeMonitor slvtimer(*timerSolve_);
666 #endif
667 
668  while ( numRHS2Solve > 0 ) {
669 
670  // Reset the active / converged vectors from this block
671  std::vector<int> convRHSIdx;
672  std::vector<int> currRHSIdx( currIdx );
673  currRHSIdx.resize(numCurrRHS);
674 
675  // Reset the number of iterations.
676  block_cg_iter->resetNumIters();
677 
678  // Reset the number of calls that the status test output knows about.
679  outputTest_->resetNumCalls();
680 
681  // Get the current residual for this block of linear systems.
682  Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
683 
684  // Get a new state struct and initialize the solver.
686  newState.R = R_0;
687  block_cg_iter->initializeCG(newState);
688 
689  while(1) {
690 
691  // tell block_gmres_iter to iterate
692  try {
693  block_cg_iter->iterate();
694 
696  //
697  // check convergence first
698  //
700  if ( convTest_->getStatus() == Passed ) {
701 
702  // Figure out which linear systems converged.
703  std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
704 
705  // If the number of converged linear systems is equal to the
706  // number of current linear systems, then we are done with this block.
707  if (convIdx.size() == currRHSIdx.size())
708  break; // break from while(1){block_cg_iter->iterate()}
709 
710  // Inform the linear problem that we are finished with this current linear system.
711  problem_->setCurrLS();
712 
713  // Reset currRHSIdx to have the right-hand sides that are left to converge for this block.
714  int have = 0;
715  std::vector<int> unconvIdx(currRHSIdx.size());
716  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
717  bool found = false;
718  for (unsigned int j=0; j<convIdx.size(); ++j) {
719  if (currRHSIdx[i] == convIdx[j]) {
720  found = true;
721  break;
722  }
723  }
724  if (!found) {
725  currIdx2[have] = currIdx2[i];
726  currRHSIdx[have++] = currRHSIdx[i];
727  }
728  }
729  currRHSIdx.resize(have);
730  currIdx2.resize(have);
731 
732  // Set the remaining indices after deflation.
733  problem_->setLSIndex( currRHSIdx );
734 
735  // Get the current residual vector.
736  std::vector<MagnitudeType> norms;
737  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
738  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
739 
740  // Set the new state and initialize the solver.
742  defstate.R = R_0;
743  block_cg_iter->initializeCG(defstate);
744  }
745 
747  //
748  // check for maximum iterations
749  //
751  else if ( maxIterTest_->getStatus() == Passed ) {
752  // we don't have convergence
753  isConverged = false;
754  break; // break from while(1){block_cg_iter->iterate()}
755  }
756 
758  //
759  // we returned from iterate(), but none of our status tests Passed.
760  // something is wrong, and it is probably our fault.
761  //
763 
764  else {
765  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
766  "Belos::PseudoBlockStochasticCGSolMgr::solve(): Invalid return from PseudoBlockStochasticCGIter::iterate().");
767  }
768  }
769  catch (const std::exception &e) {
770  printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockStochasticCGIter::iterate() at iteration "
771  << block_cg_iter->getNumIters() << std::endl
772  << e.what() << std::endl;
773  throw;
774  }
775  }
776 
777  // Inform the linear problem that we are finished with this block linear system.
778  problem_->setCurrLS();
779 
780  // Update indices for the linear systems to be solved.
781  startPtr += numCurrRHS;
782  numRHS2Solve -= numCurrRHS;
783 
784  if ( numRHS2Solve > 0 ) {
785 
786  numCurrRHS = numRHS2Solve;
787  currIdx.resize( numCurrRHS );
788  currIdx2.resize( numCurrRHS );
789  for (int i=0; i<numCurrRHS; ++i)
790  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
791 
792  // Set the next indices.
793  problem_->setLSIndex( currIdx );
794  }
795  else {
796  currIdx.resize( numRHS2Solve );
797  }
798 
799  }// while ( numRHS2Solve > 0 )
800 
801  }
802 
803  // get the final stochastic vector
804  Y_=block_cg_iter->getStochasticVector();
805 
806 
807  // print final summary
808  sTest_->print( printer_->stream(FinalSummary) );
809 
810  // print timing information
811 #ifdef BELOS_TEUCHOS_TIME_MONITOR
812  // Calling summarize() can be expensive, so don't call unless the
813  // user wants to print out timing details. summarize() will do all
814  // the work even if it's passed a "black hole" output stream.
815  if (verbosity_ & TimingDetails)
816  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
817 #endif
818 
819  // get iteration information for this solve
820  numIters_ = maxIterTest_->getNumIters();
821 
822  if (!isConverged ) {
823  return Unconverged; // return from PseudoBlockStochasticCGSolMgr::solve()
824  }
825  return Converged; // return from PseudoBlockStochasticCGSolMgr::solve()
826 }
827 
828 // This method requires the solver manager to return a std::string that describes itself.
829 template<class ScalarType, class MV, class OP>
831 {
832  std::ostringstream oss;
833  oss << "Belos::PseudoBlockStochasticCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
834  oss << "{";
835  oss << "}";
836  return oss.str();
837 }
838 
839 } // end Belos namespace
840 
841 #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
Collection of types and exceptions used within the Belos solvers.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< OutputManager< ScalarType > > printer_
Class which manages the output and verbosity of the Belos solvers.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:120
Teuchos::RCP< const Teuchos::ParameterList > validParams_
List of valid parameters and their default values.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
int getNumIters() const override
Get the iteration count for the most recent call to solve().
std::string description() const override
Method to return description of the block CG solver manager.
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:293
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
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.
Belos::StatusTest for logically combining several status tests.
A Belos::StatusTest class for specifying a maximum number of iterations.
Belos concrete class for performing the stochastic pseudo-block CG iteration.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
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.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Get a parameter list containing the valid parameters for this object.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params) override
Set the parameters the solver manager should use to solve the linear problem.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::RCP< const MV > R
The current residual.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< Teuchos::ParameterList > params_
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
PseudoBlockStochasticCGSolMgr()
Empty constructor for BlockStochasticCGSolMgr. This constructor takes no arguments and sets the defau...
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.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver&#39;s iterate() routine unti...
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
This class implements the stochastic pseudo-block CG iteration, where the basic stochastic CG algorit...
The Belos::PseudoBlockStochasticCGSolMgr provides a powerful and fully-featured solver manager over t...
Teuchos::RCP< MV > getStochasticVector()
Get a copy of the final stochastic vector.
A class for extending the status testing capabilities of Belos via logical combinations.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:283
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
PseudoBlockStochasticCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Structure to contain pointers to CGIteration state variables.