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  {}
133  BlockCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem,
134  const Teuchos::RCP<Teuchos::ParameterList>& pl) :
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:
171  typedef Teuchos::ScalarTraits<ScalarType> SCT;
172  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
173  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
174 
175  public:
176 
178 
179 
185  BlockCGSolMgr();
186 
223  BlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
224  const Teuchos::RCP<Teuchos::ParameterList> &pl );
225 
227  virtual ~BlockCGSolMgr() {};
228 
230  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
231  return Teuchos::rcp(new BlockCGSolMgr<ScalarType,MV,OP>);
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 
248  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
249 
255  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
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 
334  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
335 
337  Teuchos::RCP<OutputManager<ScalarType> > printer_;
339  Teuchos::RCP<std::ostream> outputStream_;
340 
345  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
346 
348  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
349 
351  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
352 
354  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
355 
357  Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
358 
360  Teuchos::RCP<Teuchos::ParameterList> params_;
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 
410  Teuchos::RCP<Teuchos::Time> timerSolve_;
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>
443 BlockCGSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
444  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
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
478 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
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>
765 Teuchos::RCP<const Teuchos::ParameterList>
767 {
768  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
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 
832  Teuchos::BLAS<int,ScalarType> blas;
833  Teuchos::LAPACK<int,ScalarType> lapack;
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.
872  Teuchos::ParameterList plist;
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 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
Structure to contain pointers to CGIteration state variables.
Belos concrete class for performing the conjugate-gradient (CG) iteration.
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.
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.
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...
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.
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: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 ...
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
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.
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...