Belos Package Browser (Single Doxygen Collection)  Development
 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 //
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_BLOCK_CG_SOLMGR_HPP
43 #define BELOS_BLOCK_CG_SOLMGR_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
55 #include "BelosCGIter.hpp"
56 #include "BelosCGSingleRedIter.hpp"
57 #include "BelosBlockCGIter.hpp"
61 #include "BelosStatusTestCombo.hpp"
63 #include "BelosOutputManager.hpp"
64 #include "Teuchos_BLAS.hpp"
65 #include "Teuchos_LAPACK.hpp"
66 #ifdef BELOS_TEUCHOS_TIME_MONITOR
67 # include "Teuchos_TimeMonitor.hpp"
68 #endif
69 #if defined(HAVE_TEUCHOSCORE_CXX11)
70 # include <type_traits>
71 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
72 #include <algorithm>
73 
90 namespace Belos {
91 
93 
94 
102  BlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
103  {}};
104 
111  class BlockCGSolMgrOrthoFailure : public BelosError {public:
112  BlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
113  {}};
114 
115 
116  template<class ScalarType, class MV, class OP,
117  const bool lapackSupportsScalarType =
120  public Details::SolverManagerRequiresLapack<ScalarType,MV,OP>
121  {
122  static const bool requiresLapack =
124  typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
126 
127  public:
129  base_type ()
130  {}
133  base_type ()
134  {}
135  virtual ~BlockCGSolMgr () {}
136  };
137 
138 
139  // Partial specialization for ScalarType types for which
140  // Teuchos::LAPACK has a valid implementation. This contains the
141  // actual working implementation of BlockCGSolMgr.
142  template<class ScalarType, class MV, class OP>
143  class BlockCGSolMgr<ScalarType, MV, OP, true> :
144  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
145  {
146  // This partial specialization is already chosen for those scalar types
147  // that support lapack, so we don't need to have an additional compile-time
148  // check that the scalar type is float/double/complex.
149 // #if defined(HAVE_TEUCHOSCORE_CXX11)
150 // # if defined(HAVE_TEUCHOS_COMPLEX)
151 // static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
152 // std::is_same<ScalarType, std::complex<double> >::value ||
153 // std::is_same<ScalarType, float>::value ||
154 // std::is_same<ScalarType, double>::value,
155 // "Belos::GCRODRSolMgr: ScalarType must be one of the four "
156 // "types (S,D,C,Z) supported by LAPACK.");
157 // # else
158 // static_assert (std::is_same<ScalarType, float>::value ||
159 // std::is_same<ScalarType, double>::value,
160 // "Belos::GCRODRSolMgr: ScalarType must be float or double. "
161 // "Complex arithmetic support is currently disabled. To "
162 // "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
163 // # endif // defined(HAVE_TEUCHOS_COMPLEX)
164 // #endif // defined(HAVE_TEUCHOSCORE_CXX11)
165 
166  private:
172 
173  public:
174 
176 
177 
183  BlockCGSolMgr();
184 
223 
225  virtual ~BlockCGSolMgr() {};
226 
230  }
232 
234 
235 
236  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
237  return *problem_;
238  }
239 
242  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
243 
247 
254  return Teuchos::tuple(timerSolve_);
255  }
256 
262  MagnitudeType achievedTol() const override {
263  return achievedTol_;
264  }
265 
267  int getNumIters() const override {
268  return numIters_;
269  }
270 
273  bool isLOADetected() const override { return false; }
274 
276 
278 
279 
281  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
282 
284  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
285 
287 
289 
290 
294  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
296 
298 
299 
317  ReturnType solve() override;
318 
320 
323 
325  std::string description() const override;
326 
328 
329  private:
330 
333 
338 
344 
347 
350 
353 
356 
359 
360  //
361  // Default solver parameters.
362  //
363  static constexpr int maxIters_default_ = 1000;
364  static constexpr bool adaptiveBlockSize_default_ = true;
365  static constexpr bool showMaxResNormOnly_default_ = false;
366  static constexpr bool useSingleReduction_default_ = false;
367  static constexpr int blockSize_default_ = 1;
368  static constexpr int verbosity_default_ = Belos::Errors;
369  static constexpr int outputStyle_default_ = Belos::General;
370  static constexpr int outputFreq_default_ = -1;
371  static constexpr const char * resNorm_default_ = "TwoNorm";
372  static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ = false;
373  static constexpr const char * resScale_default_ = "Norm of Initial Residual";
374  static constexpr const char * label_default_ = "Belos";
375  static constexpr const char * orthoType_default_ = "ICGS";
376  static constexpr bool assertPositiveDefiniteness_default_ = true;
377  static constexpr std::ostream * outputStream_default_ = &std::cout;
378 
379  //
380  // Current solver parameters and other values.
381  //
382 
385 
388 
395 
398 
401 
403  int blockSize_, verbosity_, outputStyle_, outputFreq_;
404  bool adaptiveBlockSize_, showMaxResNormOnly_, useSingleReduction_;
405  std::string orthoType_, resScale_;
408 
410  std::string label_;
411 
414 
416  bool isSet_;
417  };
418 
419 
420 // Empty Constructor
421 template<class ScalarType, class MV, class OP>
423  outputStream_(Teuchos::rcp(outputStream_default_,false)),
424  convtol_(DefaultSolverParameters::convTol),
425  orthoKappa_(DefaultSolverParameters::orthoKappa),
426  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
427  maxIters_(maxIters_default_),
428  numIters_(0),
429  blockSize_(blockSize_default_),
430  verbosity_(verbosity_default_),
431  outputStyle_(outputStyle_default_),
432  outputFreq_(outputFreq_default_),
433  adaptiveBlockSize_(adaptiveBlockSize_default_),
434  showMaxResNormOnly_(showMaxResNormOnly_default_),
435  useSingleReduction_(useSingleReduction_default_),
436  orthoType_(orthoType_default_),
437  resScale_(resScale_default_),
438  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
439  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
440  label_(label_default_),
441  isSet_(false)
442 {}
443 
444 
445 // Basic Constructor
446 template<class ScalarType, class MV, class OP>
450  problem_(problem),
451  outputStream_(Teuchos::rcp(outputStream_default_,false)),
452  convtol_(DefaultSolverParameters::convTol),
453  orthoKappa_(DefaultSolverParameters::orthoKappa),
454  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
455  maxIters_(maxIters_default_),
456  numIters_(0),
457  blockSize_(blockSize_default_),
458  verbosity_(verbosity_default_),
459  outputStyle_(outputStyle_default_),
460  outputFreq_(outputFreq_default_),
461  adaptiveBlockSize_(adaptiveBlockSize_default_),
462  showMaxResNormOnly_(showMaxResNormOnly_default_),
463  useSingleReduction_(useSingleReduction_default_),
464  orthoType_(orthoType_default_),
465  resScale_(resScale_default_),
466  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
467  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
468  label_(label_default_),
469  isSet_(false)
470 {
471  TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
472  "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
473 
474  // If the user passed in a nonnull parameter list, set parameters.
475  // Otherwise, the next solve() call will use default parameters,
476  // unless the user calls setParameters() first.
477  if (! pl.is_null()) {
478  setParameters (pl);
479  }
480 }
481 
482 template<class ScalarType, class MV, class OP>
483 void
486 {
487  // Create the internal parameter list if one doesn't already exist.
488  if (params_ == Teuchos::null) {
489  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
490  }
491  else {
492  params->validateParameters(*getValidParameters());
493  }
494 
495  // Check for maximum number of iterations
496  if (params->isParameter("Maximum Iterations")) {
497  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
498 
499  // Update parameter in our list and in status test.
500  params_->set("Maximum Iterations", maxIters_);
501  if (maxIterTest_!=Teuchos::null)
502  maxIterTest_->setMaxIters( maxIters_ );
503  }
504 
505  // Check for blocksize
506  if (params->isParameter("Block Size")) {
507  blockSize_ = params->get("Block Size",blockSize_default_);
508  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
509  "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
510 
511  // Update parameter in our list.
512  params_->set("Block Size", blockSize_);
513  }
514 
515  // Check if the blocksize should be adaptive
516  if (params->isParameter("Adaptive Block Size")) {
517  adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
518 
519  // Update parameter in our list.
520  params_->set("Adaptive Block Size", adaptiveBlockSize_);
521  }
522 
523  // Check if the user is requesting the single-reduction version of CG (only for blocksize == 1)
524  if (params->isParameter("Use Single Reduction")) {
525  useSingleReduction_ = params->get("Use Single Reduction", useSingleReduction_default_);
526  if (useSingleReduction_)
527  foldConvergenceDetectionIntoAllreduce_ = params->get("Fold Convergence Detection Into Allreduce",
528  foldConvergenceDetectionIntoAllreduce_default_);
529  }
530 
531  // Check to see if the timer label changed.
532  if (params->isParameter("Timer Label")) {
533  std::string tempLabel = params->get("Timer Label", label_default_);
534 
535  // Update parameter in our list and solver timer
536  if (tempLabel != label_) {
537  label_ = tempLabel;
538  params_->set("Timer Label", label_);
539  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
540 #ifdef BELOS_TEUCHOS_TIME_MONITOR
541  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
542 #endif
543  if (ortho_ != Teuchos::null) {
544  ortho_->setLabel( label_ );
545  }
546  }
547  }
548 
549  // Check for a change in verbosity level
550  if (params->isParameter("Verbosity")) {
551  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
552  verbosity_ = params->get("Verbosity", verbosity_default_);
553  } else {
554  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
555  }
556 
557  // Update parameter in our list.
558  params_->set("Verbosity", verbosity_);
559  if (printer_ != Teuchos::null)
560  printer_->setVerbosity(verbosity_);
561  }
562 
563  // Check for a change in output style
564  if (params->isParameter("Output Style")) {
565  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
566  outputStyle_ = params->get("Output Style", outputStyle_default_);
567  } else {
568  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
569  }
570 
571  // Update parameter in our list.
572  params_->set("Output Style", outputStyle_);
573  outputTest_ = Teuchos::null;
574  }
575 
576  // output stream
577  if (params->isParameter("Output Stream")) {
578  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
579 
580  // Update parameter in our list.
581  params_->set("Output Stream", outputStream_);
582  if (printer_ != Teuchos::null)
583  printer_->setOStream( outputStream_ );
584  }
585 
586  // frequency level
587  if (verbosity_ & Belos::StatusTestDetails) {
588  if (params->isParameter("Output Frequency")) {
589  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
590  }
591 
592  // Update parameter in out list and output status test.
593  params_->set("Output Frequency", outputFreq_);
594  if (outputTest_ != Teuchos::null)
595  outputTest_->setOutputFrequency( outputFreq_ );
596  }
597 
598  // Create output manager if we need to.
599  if (printer_ == Teuchos::null) {
600  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
601  }
602 
603  // Check if the orthogonalization changed.
604  bool changedOrthoType = false;
605  if (params->isParameter("Orthogonalization")) {
606  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
607  if (tempOrthoType != orthoType_) {
608  orthoType_ = tempOrthoType;
609  changedOrthoType = true;
610  }
611  }
612  params_->set("Orthogonalization", orthoType_);
613 
614  // Check which orthogonalization constant to use.
615  if (params->isParameter("Orthogonalization Constant")) {
616  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
617  orthoKappa_ = params->get ("Orthogonalization Constant",
618  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
619  }
620  else {
621  orthoKappa_ = params->get ("Orthogonalization Constant",
623  }
624 
625  // Update parameter in our list.
626  params_->set("Orthogonalization Constant",orthoKappa_);
627  if (orthoType_=="DGKS") {
628  if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
629  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
630  }
631  }
632  }
633 
634  // Create orthogonalization manager if we need to.
635  if (ortho_ == Teuchos::null || changedOrthoType) {
637  Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
638  if (orthoType_=="DGKS" && orthoKappa_ > 0) {
639  paramsOrtho->set ("depTol", orthoKappa_ );
640  }
641 
642  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
643  }
644 
645  // Convergence
646  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
647  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
648 
649  // Check for convergence tolerance
650  if (params->isParameter("Convergence Tolerance")) {
651  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
652  convtol_ = params->get ("Convergence Tolerance",
653  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
654  }
655  else {
656  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
657  }
658 
659  // Update parameter in our list and residual tests.
660  params_->set("Convergence Tolerance", convtol_);
661  if (convTest_ != Teuchos::null)
662  convTest_->setTolerance( convtol_ );
663  }
664 
665  if (params->isParameter("Show Maximum Residual Norm Only")) {
666  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
667 
668  // Update parameter in our list and residual tests
669  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
670  if (convTest_ != Teuchos::null)
671  convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
672  }
673 
674  // Check for a change in scaling, if so we need to build new residual tests.
675  bool newResTest = false;
676  {
677  std::string tempResScale = resScale_;
678  if (params->isParameter ("Implicit Residual Scaling")) {
679  tempResScale = params->get<std::string> ("Implicit Residual Scaling");
680  }
681 
682  // Only update the scaling if it's different.
683  if (resScale_ != tempResScale) {
684  const Belos::ScaleType resScaleType =
685  convertStringToScaleType (tempResScale);
686  resScale_ = tempResScale;
687 
688  // Update parameter in our list and residual tests
689  params_->set ("Implicit Residual Scaling", resScale_);
690 
691  if (! convTest_.is_null ()) {
692  try {
693  NormType normType = Belos::TwoNorm;
694  if (params->isParameter("Residual Norm")) {
695  if (params->isType<std::string> ("Residual Norm")) {
696  normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
697  }
698  }
699  convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
700  convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
701  }
702  catch (std::exception& e) {
703  // Make sure the convergence test gets constructed again.
704  newResTest = true;
705  }
706  }
707  }
708  }
709 
710  // Create status tests if we need to.
711 
712  // Basic test checks maximum iterations and native residual.
713  if (maxIterTest_ == Teuchos::null)
714  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
715 
716  // Implicit residual test, using the native residual to determine if convergence was achieved.
717  if (convTest_.is_null () || newResTest) {
718 
719  NormType normType = Belos::TwoNorm;
720  if (params->isParameter("Residual Norm")) {
721  if (params->isType<std::string> ("Residual Norm")) {
722  normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
723  }
724  }
725 
726  convTest_ = rcp (new StatusTestResNorm_t (convtol_, 1, showMaxResNormOnly_));
727  convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
728  convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
729  }
730 
731  if (sTest_.is_null () || newResTest)
732  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
733 
734  if (outputTest_.is_null () || newResTest) {
735 
736  // Create the status test output class.
737  // This class manages and formats the output from the status test.
738  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
739  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
740 
741  // Set the solver string for the output test
742  std::string solverDesc = " Block CG ";
743  outputTest_->setSolverDesc( solverDesc );
744 
745  }
746 
747  // BelosCgIter accepts a parameter specifying whether to assert for the positivity of p^H*A*p in the CG iteration
748  if (params->isParameter("Assert Positive Definiteness")) {
749  assertPositiveDefiniteness_ = Teuchos::getParameter<bool>(*params,"Assert Positive Definiteness");
750  params_->set("Assert Positive Definiteness", assertPositiveDefiniteness_);
751  }
752 
753  // Create the timer if we need to.
754  if (timerSolve_ == Teuchos::null) {
755  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
756 #ifdef BELOS_TEUCHOS_TIME_MONITOR
757  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
758 #endif
759  }
760 
761  // Inform the solver manager that the current parameters were set.
762  isSet_ = true;
763 }
764 
765 
766 template<class ScalarType, class MV, class OP>
769 {
771 
772  // Set all the valid parameters and their default values.
773  if(is_null(validPL)) {
774  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
775  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
776  "The relative residual tolerance that needs to be achieved by the\n"
777  "iterative solver in order for the linear system to be declared converged.");
778  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
779  "The maximum number of block iterations allowed for each\n"
780  "set of RHS solved.");
781  pl->set("Block Size", static_cast<int>(blockSize_default_),
782  "The number of vectors in each block.");
783  pl->set("Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
784  "Whether the solver manager should adapt to the block size\n"
785  "based on the number of RHS to solve.");
786  pl->set("Verbosity", static_cast<int>(verbosity_default_),
787  "What type(s) of solver information should be outputted\n"
788  "to the output stream.");
789  pl->set("Output Style", static_cast<int>(outputStyle_default_),
790  "What style is used for the solver information outputted\n"
791  "to the output stream.");
792  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
793  "How often convergence information should be outputted\n"
794  "to the output stream.");
795  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
796  "A reference-counted pointer to the output stream where all\n"
797  "solver output is sent.");
798  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
799  "When convergence information is printed, only show the maximum\n"
800  "relative residual norm when the block size is greater than one.");
801  pl->set("Use Single Reduction", static_cast<bool>(useSingleReduction_default_),
802  "Use single reduction iteration when the block size is one.");
803  pl->set("Implicit Residual Scaling", resScale_default_,
804  "The type of scaling used in the residual convergence test.");
805  pl->set("Timer Label", static_cast<const char *>(label_default_),
806  "The string to use as a prefix for the timer labels.");
807  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
808  "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
809  pl->set("Assert Positive Definiteness",static_cast<bool>(assertPositiveDefiniteness_default_),
810  "Assert for positivity of p^H*A*p in CG iteration.");
811  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
812  "The constant used by DGKS orthogonalization to determine\n"
813  "whether another step of classical Gram-Schmidt is necessary.");
814  pl->set("Residual Norm",static_cast<const char *>(resNorm_default_),
815  "Norm used for the convergence check on the residual.");
816  pl->set("Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
817  "Merge the allreduce for convergence detection with the one for CG.\n"
818  "This saves one all-reduce, but incurs more computation.");
819  validPL = pl;
820  }
821  return validPL;
822 }
823 
824 
825 // solve()
826 template<class ScalarType, class MV, class OP>
828  using Teuchos::RCP;
829  using Teuchos::rcp;
830  using Teuchos::rcp_const_cast;
831  using Teuchos::rcp_dynamic_cast;
832 
833  // Set the current parameters if they were not set before. NOTE:
834  // This may occur if the user generated the solver manager with the
835  // default constructor and then didn't set any parameters using
836  // setParameters().
837  if (!isSet_) {
838  setParameters(Teuchos::parameterList(*getValidParameters()));
839  }
840 
843 
844  TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(),
846  "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
847  "has not been called.");
848 
849  // Create indices for the linear systems to be solved.
850  int startPtr = 0;
851  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
852  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
853 
854  std::vector<int> currIdx, currIdx2;
855  // If an adaptive block size is allowed then only the linear
856  // systems that need to be solved are solved. Otherwise, the index
857  // set is generated that informs the linear problem that some
858  // linear systems are augmented.
859  if ( adaptiveBlockSize_ ) {
860  blockSize_ = numCurrRHS;
861  currIdx.resize( numCurrRHS );
862  currIdx2.resize( numCurrRHS );
863  for (int i=0; i<numCurrRHS; ++i)
864  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
865 
866  }
867  else {
868  currIdx.resize( blockSize_ );
869  currIdx2.resize( blockSize_ );
870  for (int i=0; i<numCurrRHS; ++i)
871  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
872  for (int i=numCurrRHS; i<blockSize_; ++i)
873  { currIdx[i] = -1; currIdx2[i] = i; }
874  }
875 
876  // Inform the linear problem of the current linear system to solve.
877  problem_->setLSIndex( currIdx );
878 
880  // Set up the parameter list for the Iteration subclass.
882  plist.set("Block Size",blockSize_);
883 
884  // Reset the output status test (controls all the other status tests).
885  outputTest_->reset();
886 
887  // Assume convergence is achieved, then let any failed convergence
888  // set this to false. "Innocent until proven guilty."
889  bool isConverged = true;
890 
892  // Set up the BlockCG Iteration subclass.
893 
894  plist.set("Assert Positive Definiteness", assertPositiveDefiniteness_);
895 
896  RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
897  if (blockSize_ == 1) {
898  // Standard (nonblock) CG is faster for the special case of a
899  // block size of 1. A single reduction iteration can also be used
900  // if collectives are more expensive than vector updates.
901  if (useSingleReduction_) {
902  plist.set("Fold Convergence Detection Into Allreduce",
903  foldConvergenceDetectionIntoAllreduce_);
904  block_cg_iter =
905  rcp (new CGSingleRedIter<ScalarType,MV,OP> (problem_, printer_,
906  outputTest_, convTest_, plist));
907  }
908  else {
909  block_cg_iter =
910  rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_,
911  outputTest_, plist));
912  }
913  } else {
914  block_cg_iter =
915  rcp (new BlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_,
916  ortho_, plist));
917  }
918 
919 
920  // Enter solve() iterations
921  {
922 #ifdef BELOS_TEUCHOS_TIME_MONITOR
923  Teuchos::TimeMonitor slvtimer(*timerSolve_);
924 #endif
925 
926  while ( numRHS2Solve > 0 ) {
927  //
928  // Reset the active / converged vectors from this block
929  std::vector<int> convRHSIdx;
930  std::vector<int> currRHSIdx( currIdx );
931  currRHSIdx.resize(numCurrRHS);
932 
933  // Reset the number of iterations.
934  block_cg_iter->resetNumIters();
935 
936  // Reset the number of calls that the status test output knows about.
937  outputTest_->resetNumCalls();
938 
939  // Get the current residual for this block of linear systems.
940  RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
941 
942  // Set the new state and initialize the solver.
944  newstate.R = R_0;
945  block_cg_iter->initializeCG(newstate);
946 
947  while(1) {
948 
949  // tell block_cg_iter to iterate
950  try {
951  block_cg_iter->iterate();
952  //
953  // Check whether any of the linear systems converged.
954  //
955  if (convTest_->getStatus() == Passed) {
956  // At least one of the linear system(s) converged.
957  //
958  // Get the column indices of the linear systems that converged.
959  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
960  std::vector<int> convIdx =
961  rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
962 
963  // If the number of converged linear systems equals the
964  // number of linear systems currently being solved, then
965  // we are done with this block.
966  if (convIdx.size() == currRHSIdx.size())
967  break; // break from while(1){block_cg_iter->iterate()}
968 
969  // Inform the linear problem that we are finished with
970  // this current linear system.
971  problem_->setCurrLS();
972 
973  // Reset currRHSIdx to contain the right-hand sides that
974  // are left to converge for this block.
975  int have = 0;
976  std::vector<int> unconvIdx(currRHSIdx.size());
977  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
978  bool found = false;
979  for (unsigned int j=0; j<convIdx.size(); ++j) {
980  if (currRHSIdx[i] == convIdx[j]) {
981  found = true;
982  break;
983  }
984  }
985  if (!found) {
986  currIdx2[have] = currIdx2[i];
987  currRHSIdx[have++] = currRHSIdx[i];
988  }
989  else {
990  }
991  }
992  currRHSIdx.resize(have);
993  currIdx2.resize(have);
994 
995  // Set the remaining indices after deflation.
996  problem_->setLSIndex( currRHSIdx );
997 
998  // Get the current residual vector.
999  std::vector<MagnitudeType> norms;
1000  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
1001  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
1002 
1003  // Set the new blocksize for the solver.
1004  block_cg_iter->setBlockSize( have );
1005 
1006  // Set the new state and initialize the solver.
1008  defstate.R = R_0;
1009  block_cg_iter->initializeCG(defstate);
1010  }
1011  //
1012  // None of the linear systems converged. Check whether the
1013  // maximum iteration count was reached.
1014  //
1015  else if (maxIterTest_->getStatus() == Passed) {
1016  isConverged = false; // None of the linear systems converged.
1017  break; // break from while(1){block_cg_iter->iterate()}
1018  }
1019  //
1020  // iterate() returned, but none of our status tests Passed.
1021  // This indicates a bug.
1022  //
1023  else {
1024  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1025  "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
1026  "the maximum iteration count test passed. Please report this bug "
1027  "to the Belos developers.");
1028  }
1029  }
1030  catch (const std::exception &e) {
1031  std::ostream& err = printer_->stream (Errors);
1032  err << "Error! Caught std::exception in CGIteration::iterate() at "
1033  << "iteration " << block_cg_iter->getNumIters() << std::endl
1034  << e.what() << std::endl;
1035  throw;
1036  }
1037  }
1038 
1039  // Inform the linear problem that we are finished with this
1040  // block linear system.
1041  problem_->setCurrLS();
1042 
1043  // Update indices for the linear systems to be solved.
1044  startPtr += numCurrRHS;
1045  numRHS2Solve -= numCurrRHS;
1046  if ( numRHS2Solve > 0 ) {
1047  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1048 
1049  if ( adaptiveBlockSize_ ) {
1050  blockSize_ = numCurrRHS;
1051  currIdx.resize( numCurrRHS );
1052  currIdx2.resize( numCurrRHS );
1053  for (int i=0; i<numCurrRHS; ++i)
1054  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1055  }
1056  else {
1057  currIdx.resize( blockSize_ );
1058  currIdx2.resize( blockSize_ );
1059  for (int i=0; i<numCurrRHS; ++i)
1060  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1061  for (int i=numCurrRHS; i<blockSize_; ++i)
1062  { currIdx[i] = -1; currIdx2[i] = i; }
1063  }
1064  // Set the next indices.
1065  problem_->setLSIndex( currIdx );
1066 
1067  // Set the new blocksize for the solver.
1068  block_cg_iter->setBlockSize( blockSize_ );
1069  }
1070  else {
1071  currIdx.resize( numRHS2Solve );
1072  }
1073 
1074  }// while ( numRHS2Solve > 0 )
1075 
1076  }
1077 
1078  // print final summary
1079  sTest_->print( printer_->stream(FinalSummary) );
1080 
1081  // print timing information
1082 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1083  // Calling summarize() requires communication in general, so don't
1084  // call it unless the user wants to print out timing details.
1085  // summarize() will do all the work even if it's passed a "black
1086  // hole" output stream.
1087  if (verbosity_ & TimingDetails) {
1088  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1089  }
1090 #endif
1091 
1092  // Save the iteration count for this solve.
1093  numIters_ = maxIterTest_->getNumIters();
1094 
1095  // Save the convergence test value ("achieved tolerance") for this solve.
1096  {
1097  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1098  // testValues is nonnull and not persistent.
1099  const std::vector<MagnitudeType>* pTestValues =
1100  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1101 
1102  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1103  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1104  "method returned NULL. Please report this bug to the Belos developers.");
1105 
1106  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1107  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1108  "method returned a vector of length zero. Please report this bug to the "
1109  "Belos developers.");
1110 
1111  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1112  // achieved tolerances for all vectors in the current solve(), or
1113  // just for the vectors from the last deflation?
1114  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1115  }
1116 
1117  if (!isConverged) {
1118  return Unconverged; // return from BlockCGSolMgr::solve()
1119  }
1120  return Converged; // return from BlockCGSolMgr::solve()
1121 }
1122 
1123 // This method requires the solver manager to return a std::string that describes itself.
1124 template<class ScalarType, class MV, class OP>
1126 {
1127  std::ostringstream oss;
1128  oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1129  oss << "{";
1130  oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
1131  oss << "}";
1132  return oss.str();
1133 }
1134 
1135 } // end Belos namespace
1136 
1137 #endif /* BELOS_BLOCK_CG_SOLMGR_HPP */
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:299
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
Teuchos::RCP< const MV > R
The current residual.
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...
int maxIters_
Maximum iteration count (read from parameter list).
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
bool isSet_
Whether or not the parameters have been set (via setParameters()).
Teuchos::RCP< Teuchos::ParameterList > params_
Current parameter list.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Output &quot;status test&quot; that controls all the other status tests.
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:120
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:78
T & get(ParameterList &l, const std::string &name)
Structure to contain pointers to CGIteration state variables.
Belos concrete class for performing the conjugate-gradient (CG) iteration.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
#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.
std::string label_
Prefix label for all the timers.
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:293
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.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
BlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Maximum iteration count stopping criterion.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Aggregate stopping criterion.
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, requiresLapack > base_type
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:206
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)
bool is_null(const RCP< T > &p)
NormType convertStringToNormType(const std::string &normType)
Convert the given string to its NormType enum value.
Definition: BelosTypes.cpp:88
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
BlockCGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonor...
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...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Convergence stopping criterion.
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The linear problem to solve.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
BlockCGSolMgrOrthoFailure(const std::string &what_arg)
NormType
The type of vector norm to compute.
Definition: BelosTypes.hpp:97
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.
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
bool isType(const std::string &name) const
Teuchos::RCP< Teuchos::Time > timerSolve_
Solve timer.
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:60
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:283
This class implements the preconditioned single-reduction Conjugate Gradient (CG) iteration...
Teuchos::RCP< OutputManager< ScalarType > > printer_
Output manager, that handles printing of different kinds of messages.
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...
int numIters_
Number of iterations taken by the last solve() invocation.
Teuchos::ScalarTraits< MagnitudeType > MT
static const bool requiresLapack
Teuchos::RCP< std::ostream > outputStream_
Output stream to which the output manager prints.
MagnitudeType orthoKappa_
Orthogonalization parameter (read from parameter list).
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...
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.