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