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

Generated for Belos by doxygen 1.8.14