Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosBlockCGSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_BLOCK_CG_SOLMGR_HPP
11 #define BELOS_BLOCK_CG_SOLMGR_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 
20 #include "BelosLinearProblem.hpp"
21 #include "BelosSolverManager.hpp"
22 
23 #include "BelosCGIter.hpp"
24 #include "BelosCGSingleRedIter.hpp"
25 #include "BelosBlockCGIter.hpp"
29 #include "BelosStatusTestCombo.hpp"
31 #include "BelosOutputManager.hpp"
32 #include "Teuchos_LAPACK.hpp"
33 #ifdef BELOS_TEUCHOS_TIME_MONITOR
34 # include "Teuchos_TimeMonitor.hpp"
35 #endif
36 #if defined(HAVE_TEUCHOSCORE_CXX11)
37 # include <type_traits>
38 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
39 #include <algorithm>
40 
60 namespace Belos {
61 
63 
64 
72  BlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
73  {}};
74 
75  template<class ScalarType, class MV, class OP,
76  const bool lapackSupportsScalarType =
78  class BlockCGSolMgr :
79  public Details::SolverManagerRequiresLapack<ScalarType,MV,OP>
80  {
81  static const bool requiresLapack =
83  typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
84  requiresLapack> base_type;
85 
86  public:
88  base_type ()
89  {}
92  base_type ()
93  {}
94  virtual ~BlockCGSolMgr () {}
95  };
96 
97 
98  // Partial specialization for ScalarType types for which
99  // Teuchos::LAPACK has a valid implementation. This contains the
100  // actual working implementation of BlockCGSolMgr.
101  template<class ScalarType, class MV, class OP>
102  class BlockCGSolMgr<ScalarType, MV, OP, true> :
103  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
104  {
105  // This partial specialization is already chosen for those scalar types
106  // that support lapack, so we don't need to have an additional compile-time
107  // check that the scalar type is float/double/complex.
108 // #if defined(HAVE_TEUCHOSCORE_CXX11)
109 // # if defined(HAVE_TEUCHOS_COMPLEX)
110 // static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
111 // std::is_same<ScalarType, std::complex<double> >::value ||
112 // std::is_same<ScalarType, float>::value ||
113 // std::is_same<ScalarType, double>::value,
114 // "Belos::GCRODRSolMgr: ScalarType must be one of the four "
115 // "types (S,D,C,Z) supported by LAPACK.");
116 // # else
117 // static_assert (std::is_same<ScalarType, float>::value ||
118 // std::is_same<ScalarType, double>::value,
119 // "Belos::GCRODRSolMgr: ScalarType must be float or double. "
120 // "Complex arithmetic support is currently disabled. To "
121 // "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
122 // # endif // defined(HAVE_TEUCHOS_COMPLEX)
123 // #endif // defined(HAVE_TEUCHOSCORE_CXX11)
124 
125  private:
129  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
131 
132  public:
133 
135 
136 
142  BlockCGSolMgr();
143 
182 
184  virtual ~BlockCGSolMgr() {};
185 
189  }
191 
193 
194 
195  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
196  return *problem_;
197  }
198 
201  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
202 
206 
213  return Teuchos::tuple(timerSolve_);
214  }
215 
221  MagnitudeType achievedTol() const override {
222  return achievedTol_;
223  }
224 
226  int getNumIters() const override {
227  return numIters_;
228  }
229 
232  bool isLOADetected() const override { return false; }
233 
235 
237 
238 
240  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
241 
243  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
244 
246 
248 
249 
253  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
255 
257 
258 
276  ReturnType solve() override;
277 
279 
282 
284  std::string description() const override;
285 
287 
288  private:
289 
292 
296  Teuchos::RCP<std::ostream> outputStream_;
297 
303 
306 
309 
312 
315 
318 
319  //
320  // Default solver parameters.
321  //
322  static constexpr int maxIters_default_ = 1000;
323  static constexpr bool adaptiveBlockSize_default_ = true;
324  static constexpr bool showMaxResNormOnly_default_ = false;
325  static constexpr bool useSingleReduction_default_ = false;
326  static constexpr int blockSize_default_ = 1;
327  static constexpr int verbosity_default_ = Belos::Errors;
328  static constexpr int outputStyle_default_ = Belos::General;
329  static constexpr int outputFreq_default_ = -1;
330  static constexpr const char * resNorm_default_ = "TwoNorm";
331  static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ = false;
332  static constexpr const char * resScale_default_ = "Norm of Initial Residual";
333  static constexpr const char * label_default_ = "Belos";
334  static constexpr const char * orthoType_default_ = "ICGS";
335  static constexpr bool assertPositiveDefiniteness_default_ = true;
336 
337  //
338  // Current solver parameters and other values.
339  //
340 
342  MagnitudeType convtol_;
343 
345  MagnitudeType orthoKappa_;
346 
352  MagnitudeType achievedTol_;
353 
355  int maxIters_;
356 
358  int numIters_;
359 
361  int blockSize_, verbosity_, outputStyle_, outputFreq_;
362  bool adaptiveBlockSize_, showMaxResNormOnly_, useSingleReduction_;
363  std::string orthoType_, resScale_;
364  bool assertPositiveDefiniteness_;
365  bool foldConvergenceDetectionIntoAllreduce_;
366 
368  std::string label_;
369 
371  Teuchos::RCP<Teuchos::Time> timerSolve_;
372 
374  bool isSet_;
375  };
376 
377 
378 // Empty Constructor
379 template<class ScalarType, class MV, class OP>
381  outputStream_(Teuchos::rcpFromRef(std::cout)),
382  convtol_(DefaultSolverParameters::convTol),
383  orthoKappa_(DefaultSolverParameters::orthoKappa),
384  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
385  maxIters_(maxIters_default_),
386  numIters_(0),
387  blockSize_(blockSize_default_),
388  verbosity_(verbosity_default_),
389  outputStyle_(outputStyle_default_),
390  outputFreq_(outputFreq_default_),
391  adaptiveBlockSize_(adaptiveBlockSize_default_),
392  showMaxResNormOnly_(showMaxResNormOnly_default_),
393  useSingleReduction_(useSingleReduction_default_),
394  orthoType_(orthoType_default_),
395  resScale_(resScale_default_),
396  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
397  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
398  label_(label_default_),
399  isSet_(false)
400 {}
401 
402 
403 // Basic Constructor
404 template<class ScalarType, class MV, class OP>
408  problem_(problem),
409  outputStream_(Teuchos::rcpFromRef(std::cout)),
410  convtol_(DefaultSolverParameters::convTol),
411  orthoKappa_(DefaultSolverParameters::orthoKappa),
412  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
413  maxIters_(maxIters_default_),
414  numIters_(0),
415  blockSize_(blockSize_default_),
416  verbosity_(verbosity_default_),
417  outputStyle_(outputStyle_default_),
418  outputFreq_(outputFreq_default_),
419  adaptiveBlockSize_(adaptiveBlockSize_default_),
420  showMaxResNormOnly_(showMaxResNormOnly_default_),
421  useSingleReduction_(useSingleReduction_default_),
422  orthoType_(orthoType_default_),
423  resScale_(resScale_default_),
424  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
425  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
426  label_(label_default_),
427  isSet_(false)
428 {
429  TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
430  "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
431 
432  // If the user passed in a nonnull parameter list, set parameters.
433  // Otherwise, the next solve() call will use default parameters,
434  // unless the user calls setParameters() first.
435  if (! pl.is_null()) {
436  setParameters (pl);
437  }
438 }
439 
440 template<class ScalarType, class MV, class OP>
441 void
444 {
445  // Create the internal parameter list if one doesn't already exist.
446  if (params_ == Teuchos::null) {
447  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
448  }
449  else {
450  params->validateParameters(*getValidParameters());
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_!=Teuchos::null)
460  maxIterTest_->setMaxIters( maxIters_ );
461  }
462 
463  // Check for blocksize
464  if (params->isParameter("Block Size")) {
465  blockSize_ = params->get("Block Size",blockSize_default_);
466  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
467  "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
468 
469  // Update parameter in our list.
470  params_->set("Block Size", blockSize_);
471  }
472 
473  // Check if the blocksize should be adaptive
474  if (params->isParameter("Adaptive Block Size")) {
475  adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
476 
477  // Update parameter in our list.
478  params_->set("Adaptive Block Size", adaptiveBlockSize_);
479  }
480 
481  // Check if the user is requesting the single-reduction version of CG (only for blocksize == 1)
482  if (params->isParameter("Use Single Reduction")) {
483  useSingleReduction_ = params->get("Use Single Reduction", useSingleReduction_default_);
484  }
485 
486  if (params->isParameter("Fold Convergence Detection Into Allreduce")) {
487  foldConvergenceDetectionIntoAllreduce_ = params->get("Fold Convergence Detection Into Allreduce",
488  foldConvergenceDetectionIntoAllreduce_default_);
489  }
490 
491  // Check to see if the timer label changed.
492  if (params->isParameter("Timer Label")) {
493  std::string tempLabel = params->get("Timer Label", label_default_);
494 
495  // Update parameter in our list and solver timer
496  if (tempLabel != label_) {
497  label_ = tempLabel;
498  params_->set("Timer Label", label_);
499  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
500 #ifdef BELOS_TEUCHOS_TIME_MONITOR
501  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
502 #endif
503  if (ortho_ != Teuchos::null) {
504  ortho_->setLabel( label_ );
505  }
506  }
507  }
508 
509  // Check for a change in verbosity level
510  if (params->isParameter("Verbosity")) {
511  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
512  verbosity_ = params->get("Verbosity", verbosity_default_);
513  } else {
514  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
515  }
516 
517  // Update parameter in our list.
518  params_->set("Verbosity", verbosity_);
519  if (printer_ != Teuchos::null)
520  printer_->setVerbosity(verbosity_);
521  }
522 
523  // Check for a change in output style
524  if (params->isParameter("Output Style")) {
525  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
526  outputStyle_ = params->get("Output Style", outputStyle_default_);
527  } else {
528  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
529  }
530 
531  // Update parameter in our list.
532  params_->set("Output Style", outputStyle_);
533  outputTest_ = Teuchos::null;
534  }
535 
536  // output stream
537  if (params->isParameter("Output Stream")) {
538  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
539 
540  // Update parameter in our list.
541  params_->set("Output Stream", outputStream_);
542  if (printer_ != Teuchos::null)
543  printer_->setOStream( outputStream_ );
544  }
545 
546  // frequency level
547  if (verbosity_ & Belos::StatusTestDetails) {
548  if (params->isParameter("Output Frequency")) {
549  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
550  }
551 
552  // Update parameter in out list and output status test.
553  params_->set("Output Frequency", outputFreq_);
554  if (outputTest_ != Teuchos::null)
555  outputTest_->setOutputFrequency( outputFreq_ );
556  }
557 
558  // Create output manager if we need to.
559  if (printer_ == Teuchos::null) {
560  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
561  }
562 
563  // Check if the orthogonalization changed.
564  bool changedOrthoType = false;
565  if (params->isParameter("Orthogonalization")) {
566  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
567  if (tempOrthoType != orthoType_) {
568  orthoType_ = tempOrthoType;
569  changedOrthoType = true;
570  }
571  }
572  params_->set("Orthogonalization", orthoType_);
573 
574  // Check which orthogonalization constant to use.
575  if (params->isParameter("Orthogonalization Constant")) {
576  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
577  orthoKappa_ = params->get ("Orthogonalization Constant",
578  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
579  }
580  else {
581  orthoKappa_ = params->get ("Orthogonalization Constant",
583  }
584 
585  // Update parameter in our list.
586  params_->set("Orthogonalization Constant",orthoKappa_);
587  if (orthoType_=="DGKS") {
588  if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
589  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
590  }
591  }
592  }
593 
594  // Create orthogonalization manager if we need to.
595  if (ortho_ == Teuchos::null || changedOrthoType) {
598  if (orthoType_=="DGKS" && orthoKappa_ > 0) {
599  paramsOrtho = Teuchos::rcp(new Teuchos::ParameterList());
600  paramsOrtho->set ("depTol", orthoKappa_ );
601  }
602 
603  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
604  }
605 
606  // Convergence
607  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
608  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
609 
610  // Check for convergence tolerance
611  if (params->isParameter("Convergence Tolerance")) {
612  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
613  convtol_ = params->get ("Convergence Tolerance",
614  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
615  }
616  else {
617  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
618  }
619 
620  // Update parameter in our list and residual tests.
621  params_->set("Convergence Tolerance", convtol_);
622  if (convTest_ != Teuchos::null)
623  convTest_->setTolerance( convtol_ );
624  }
625 
626  if (params->isParameter("Show Maximum Residual Norm Only")) {
627  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
628 
629  // Update parameter in our list and residual tests
630  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
631  if (convTest_ != Teuchos::null)
632  convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
633  }
634 
635  // Check for a change in scaling, if so we need to build new residual tests.
636  bool newResTest = false;
637  {
638  std::string tempResScale = resScale_;
639  if (params->isParameter ("Implicit Residual Scaling")) {
640  tempResScale = params->get<std::string> ("Implicit Residual Scaling");
641  }
642 
643  // Only update the scaling if it's different.
644  if (resScale_ != tempResScale) {
645  const Belos::ScaleType resScaleType =
646  convertStringToScaleType (tempResScale);
647  resScale_ = tempResScale;
648 
649  // Update parameter in our list and residual tests
650  params_->set ("Implicit Residual Scaling", resScale_);
651 
652  if (! convTest_.is_null ()) {
653  try {
654  NormType normType = Belos::TwoNorm;
655  if (params->isParameter("Residual Norm")) {
656  if (params->isType<std::string> ("Residual Norm")) {
657  normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
658  }
659  }
660  convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
661  convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
662  }
663  catch (std::exception& e) {
664  // Make sure the convergence test gets constructed again.
665  newResTest = true;
666  }
667  }
668  }
669  }
670 
671  // Create status tests if we need to.
672 
673  // Basic test checks maximum iterations and native residual.
674  if (maxIterTest_ == Teuchos::null)
675  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
676 
677  // Implicit residual test, using the native residual to determine if convergence was achieved.
678  if (convTest_.is_null () || newResTest) {
679 
680  NormType normType = Belos::TwoNorm;
681  if (params->isParameter("Residual Norm")) {
682  if (params->isType<std::string> ("Residual Norm")) {
683  normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
684  }
685  }
686 
687  convTest_ = rcp (new StatusTestResNorm_t (convtol_, 1, showMaxResNormOnly_));
688  convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
689  convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
690  }
691 
692  if (sTest_.is_null () || newResTest)
693  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
694 
695  if (outputTest_.is_null () || newResTest) {
696 
697  // Create the status test output class.
698  // This class manages and formats the output from the status test.
699  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
700  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
701 
702  // Set the solver string for the output test
703  std::string solverDesc = " Block CG ";
704  outputTest_->setSolverDesc( solverDesc );
705 
706  }
707 
708  // BelosCgIter accepts a parameter specifying whether to assert for the positivity of p^H*A*p in the CG iteration
709  if (params->isParameter("Assert Positive Definiteness")) {
710  assertPositiveDefiniteness_ = Teuchos::getParameter<bool>(*params,"Assert Positive Definiteness");
711  params_->set("Assert Positive Definiteness", assertPositiveDefiniteness_);
712  }
713 
714  // Create the timer if we need to.
715  if (timerSolve_ == Teuchos::null) {
716  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
717 #ifdef BELOS_TEUCHOS_TIME_MONITOR
718  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
719 #endif
720  }
721 
722  // Inform the solver manager that the current parameters were set.
723  isSet_ = true;
724 }
725 
726 
727 template<class ScalarType, class MV, class OP>
730 {
732 
733  // Set all the valid parameters and their default values.
734  if(is_null(validPL)) {
735  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
736  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
737  "The relative residual tolerance that needs to be achieved by the\n"
738  "iterative solver in order for the linear system to be declared converged.");
739  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
740  "The maximum number of block iterations allowed for each\n"
741  "set of RHS solved.");
742  pl->set("Block Size", static_cast<int>(blockSize_default_),
743  "The number of vectors in each block.");
744  pl->set("Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
745  "Whether the solver manager should adapt to the block size\n"
746  "based on the number of RHS to solve.");
747  pl->set("Verbosity", static_cast<int>(verbosity_default_),
748  "What type(s) of solver information should be outputted\n"
749  "to the output stream.");
750  pl->set("Output Style", static_cast<int>(outputStyle_default_),
751  "What style is used for the solver information outputted\n"
752  "to the output stream.");
753  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
754  "How often convergence information should be outputted\n"
755  "to the output stream.");
756  pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
757  "A reference-counted pointer to the output stream where all\n"
758  "solver output is sent.");
759  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
760  "When convergence information is printed, only show the maximum\n"
761  "relative residual norm when the block size is greater than one.");
762  pl->set("Use Single Reduction", static_cast<bool>(useSingleReduction_default_),
763  "Use single reduction iteration when the block size is one.");
764  pl->set("Implicit Residual Scaling", resScale_default_,
765  "The type of scaling used in the residual convergence test.");
766  pl->set("Timer Label", static_cast<const char *>(label_default_),
767  "The string to use as a prefix for the timer labels.");
768  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
769  "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
770  pl->set("Assert Positive Definiteness",static_cast<bool>(assertPositiveDefiniteness_default_),
771  "Assert for positivity of p^H*A*p in CG iteration.");
772  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
773  "The constant used by DGKS orthogonalization to determine\n"
774  "whether another step of classical Gram-Schmidt is necessary.");
775  pl->set("Residual Norm",static_cast<const char *>(resNorm_default_),
776  "Norm used for the convergence check on the residual.");
777  pl->set("Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
778  "Merge the allreduce for convergence detection with the one for CG.\n"
779  "This saves one all-reduce, but incurs more computation.");
780  validPL = pl;
781  }
782  return validPL;
783 }
784 
785 
786 // solve()
787 template<class ScalarType, class MV, class OP>
789  using Teuchos::RCP;
790  using Teuchos::rcp;
791  using Teuchos::rcp_const_cast;
792  using Teuchos::rcp_dynamic_cast;
793 
794  // Set the current parameters if they were not set before. NOTE:
795  // This may occur if the user generated the solver manager with the
796  // default constructor and then didn't set any parameters using
797  // setParameters().
798  if (!isSet_) {
799  setParameters(Teuchos::parameterList(*getValidParameters()));
800  }
801 
803 
804  TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(),
806  "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
807  "has not been called.");
808 
809  // Create indices for the linear systems to be solved.
810  int startPtr = 0;
811  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
812  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
813 
814  std::vector<int> currIdx, currIdx2;
815  // If an adaptive block size is allowed then only the linear
816  // systems that need to be solved are solved. Otherwise, the index
817  // set is generated that informs the linear problem that some
818  // linear systems are augmented.
819  if ( adaptiveBlockSize_ ) {
820  blockSize_ = numCurrRHS;
821  currIdx.resize( numCurrRHS );
822  currIdx2.resize( numCurrRHS );
823  for (int i=0; i<numCurrRHS; ++i)
824  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
825 
826  }
827  else {
828  currIdx.resize( blockSize_ );
829  currIdx2.resize( blockSize_ );
830  for (int i=0; i<numCurrRHS; ++i)
831  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
832  for (int i=numCurrRHS; i<blockSize_; ++i)
833  { currIdx[i] = -1; currIdx2[i] = i; }
834  }
835 
836  // Inform the linear problem of the current linear system to solve.
837  problem_->setLSIndex( currIdx );
838 
840  // Set up the parameter list for the Iteration subclass.
842  plist.set("Block Size",blockSize_);
843 
844  // Reset the output status test (controls all the other status tests).
845  outputTest_->reset();
846 
847  // Assume convergence is achieved, then let any failed convergence
848  // set this to false. "Innocent until proven guilty."
849  bool isConverged = true;
850 
852  // Set up the BlockCG Iteration subclass.
853 
854  plist.set("Assert Positive Definiteness", assertPositiveDefiniteness_);
855 
856  RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
857  if (blockSize_ == 1) {
858  // Standard (nonblock) CG is faster for the special case of a
859  // block size of 1. A single reduction iteration can also be used
860  // if collectives are more expensive than vector updates.
861  plist.set("Fold Convergence Detection Into Allreduce",
862  foldConvergenceDetectionIntoAllreduce_);
863  if (useSingleReduction_) {
864  block_cg_iter =
865  rcp (new CGSingleRedIter<ScalarType,MV,OP> (problem_, printer_,
866  outputTest_, convTest_, plist));
867  }
868  else {
869  block_cg_iter =
870  rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_,
871  outputTest_, convTest_, plist));
872  }
873  } else {
874  block_cg_iter =
875  rcp (new BlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_,
876  ortho_, plist));
877  }
878 
879 
880  // Enter solve() iterations
881  {
882 #ifdef BELOS_TEUCHOS_TIME_MONITOR
883  Teuchos::TimeMonitor slvtimer(*timerSolve_);
884 #endif
885 
886  while ( numRHS2Solve > 0 ) {
887  //
888  // Reset the active / converged vectors from this block
889  std::vector<int> convRHSIdx;
890  std::vector<int> currRHSIdx( currIdx );
891  currRHSIdx.resize(numCurrRHS);
892 
893  // Reset the number of iterations.
894  block_cg_iter->resetNumIters();
895 
896  // Reset the number of calls that the status test output knows about.
897  outputTest_->resetNumCalls();
898 
899  // Get the current residual for this block of linear systems.
900  RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
901 
902  // Set the new state and initialize the solver.
904  newstate.R = R_0;
905  block_cg_iter->initializeCG(newstate);
906 
907  while(1) {
908 
909  // tell block_cg_iter to iterate
910  try {
911  block_cg_iter->iterate();
912  //
913  // Check whether any of the linear systems converged.
914  //
915  if (convTest_->getStatus() == Passed) {
916  // At least one of the linear system(s) converged.
917  //
918  // Get the column indices of the linear systems that converged.
919  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
920  std::vector<int> convIdx =
921  rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
922 
923  // If the number of converged linear systems equals the
924  // number of linear systems currently being solved, then
925  // we are done with this block.
926  if (convIdx.size() == currRHSIdx.size())
927  break; // break from while(1){block_cg_iter->iterate()}
928 
929  // Inform the linear problem that we are finished with
930  // this current linear system.
931  problem_->setCurrLS();
932 
933  // Reset currRHSIdx to contain the right-hand sides that
934  // are left to converge for this block.
935  int have = 0;
936  std::vector<int> unconvIdx(currRHSIdx.size());
937  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
938  bool found = false;
939  for (unsigned int j=0; j<convIdx.size(); ++j) {
940  if (currRHSIdx[i] == convIdx[j]) {
941  found = true;
942  break;
943  }
944  }
945  if (!found) {
946  currIdx2[have] = currIdx2[i];
947  currRHSIdx[have++] = currRHSIdx[i];
948  }
949  else {
950  }
951  }
952  currRHSIdx.resize(have);
953  currIdx2.resize(have);
954 
955  // Set the remaining indices after deflation.
956  problem_->setLSIndex( currRHSIdx );
957 
958  // Get the current residual vector.
959  std::vector<MagnitudeType> norms;
960  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
961  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
962 
963  // Set the new blocksize for the solver.
964  block_cg_iter->setBlockSize( have );
965 
966  // Set the new state and initialize the solver.
968  defstate.R = R_0;
969  block_cg_iter->initializeCG(defstate);
970  }
971  //
972  // None of the linear systems converged. Check whether the
973  // maximum iteration count was reached.
974  //
975  else if (maxIterTest_->getStatus() == Passed) {
976  isConverged = false; // None of the linear systems converged.
977  break; // break from while(1){block_cg_iter->iterate()}
978  }
979  //
980  // iterate() returned, but none of our status tests Passed.
981  // This indicates a bug.
982  //
983  else {
984  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
985  "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
986  "the maximum iteration count test passed. Please report this bug "
987  "to the Belos developers.");
988  }
989  }
990  catch (const StatusTestNaNError& e) {
991  // A NaN was detected in the solver. Set the solution to zero and return unconverged.
992  achievedTol_ = MT::one();
993  Teuchos::RCP<MV> X = problem_->getLHS();
994  MVT::MvInit( *X, SCT::zero() );
995  printer_->stream(Warnings) << "Belos::BlockCGSolMgr::solve(): Warning! NaN has been detected!"
996  << std::endl;
997  return Unconverged;
998  }
999  catch (const std::exception &e) {
1000  std::ostream& err = printer_->stream (Errors);
1001  err << "Error! Caught std::exception in CGIteration::iterate() at "
1002  << "iteration " << block_cg_iter->getNumIters() << std::endl
1003  << e.what() << std::endl;
1004  throw;
1005  }
1006  }
1007 
1008  // Inform the linear problem that we are finished with this
1009  // block linear system.
1010  problem_->setCurrLS();
1011 
1012  // Update indices for the linear systems to be solved.
1013  startPtr += numCurrRHS;
1014  numRHS2Solve -= numCurrRHS;
1015  if ( numRHS2Solve > 0 ) {
1016  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1017 
1018  if ( adaptiveBlockSize_ ) {
1019  blockSize_ = numCurrRHS;
1020  currIdx.resize( numCurrRHS );
1021  currIdx2.resize( numCurrRHS );
1022  for (int i=0; i<numCurrRHS; ++i)
1023  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1024  }
1025  else {
1026  currIdx.resize( blockSize_ );
1027  currIdx2.resize( blockSize_ );
1028  for (int i=0; i<numCurrRHS; ++i)
1029  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1030  for (int i=numCurrRHS; i<blockSize_; ++i)
1031  { currIdx[i] = -1; currIdx2[i] = i; }
1032  }
1033  // Set the next indices.
1034  problem_->setLSIndex( currIdx );
1035 
1036  // Set the new blocksize for the solver.
1037  block_cg_iter->setBlockSize( blockSize_ );
1038  }
1039  else {
1040  currIdx.resize( numRHS2Solve );
1041  }
1042 
1043  }// while ( numRHS2Solve > 0 )
1044 
1045  }
1046 
1047  // print final summary
1048  sTest_->print( printer_->stream(FinalSummary) );
1049 
1050  // print timing information
1051 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1052  // Calling summarize() requires communication in general, so don't
1053  // call it unless the user wants to print out timing details.
1054  // summarize() will do all the work even if it's passed a "black
1055  // hole" output stream.
1056  if (verbosity_ & TimingDetails) {
1057  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1058  }
1059 #endif
1060 
1061  // Save the iteration count for this solve.
1062  numIters_ = maxIterTest_->getNumIters();
1063 
1064  // Save the convergence test value ("achieved tolerance") for this solve.
1065  {
1066  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1067  // testValues is nonnull and not persistent.
1068  const std::vector<MagnitudeType>* pTestValues =
1069  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1070 
1071  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1072  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1073  "method returned NULL. Please report this bug to the Belos developers.");
1074 
1075  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1076  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1077  "method returned a vector of length zero. Please report this bug to the "
1078  "Belos developers.");
1079 
1080  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1081  // achieved tolerances for all vectors in the current solve(), or
1082  // just for the vectors from the last deflation?
1083  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1084  }
1085 
1086  if (!isConverged) {
1087  return Unconverged; // return from BlockCGSolMgr::solve()
1088  }
1089  return Converged; // return from BlockCGSolMgr::solve()
1090 }
1091 
1092 // This method requires the solver manager to return a std::string that describes itself.
1093 template<class ScalarType, class MV, class OP>
1095 {
1096  std::ostringstream oss;
1097  oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1098  oss << "{";
1099  oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
1100  oss << "}";
1101  return oss.str();
1102 }
1103 
1104 } // end Belos namespace
1105 
1106 #endif /* BELOS_BLOCK_CG_SOLMGR_HPP */
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:267
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:74
Teuchos::RCP< const MV > R
The current residual.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Collection of types and exceptions used within the Belos solvers.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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.
bool is_null(const boost::shared_ptr< T > &p)
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:88
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
This class implements the preconditioned Conjugate Gradient (CG) iteration.
Definition: BelosCGIter.hpp:47
Structure to contain pointers to CGIteration state variables.
T & get(const std::string &name, T def_value)
Belos concrete class for performing the conjugate-gradient (CG) iteration.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
Belos concrete class for performing the block conjugate-gradient (CG) iteration.
An implementation of StatusTestResNorm using a family of residual norms.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
The Belos::BlockCGSolMgr provides a powerful and fully-featured solver manager over the CG and BlockC...
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:261
Belos::StatusTest class for specifying a maximum number of iterations.
static std::string name()
BlockCGSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
BlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
A factory class for generating StatusTestOutput objects.
BlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
bool isParameter(const std::string &name) const
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:174
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
Belos concrete class for performing a single-reduction conjugate-gradient (CG) iteration.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
NormType convertStringToNormType(const std::string &normType)
Convert the given string to its NormType enum value.
Definition: BelosTypes.cpp:56
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:123
int getNumIters() const override
Get the iteration count for the most recent call to solve().
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
NormType
The type of vector norm to compute.
Definition: BelosTypes.hpp:65
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.
bool isType(const std::string &name) const
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:251
This class implements the preconditioned single-reduction Conjugate Gradient (CG) iteration...
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...
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...
bool is_null() const
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
Return an instance of the specified MatOrthoManager subclass.

Generated on Mon Jul 15 2024 09:24:23 for Belos by doxygen 1.8.5