Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosBlockCGSolMgr.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
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,
127  requiresLapack> base_type;
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:
172  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
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 
339  Teuchos::RCP<std::ostream> outputStream_;
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 * resNorm_default_ = "TwoNorm";
374  static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ = false;
375  static constexpr const char * resScale_default_ = "Norm of Initial Residual";
376  static constexpr const char * label_default_ = "Belos";
377  static constexpr const char * orthoType_default_ = "DGKS";
378  static constexpr std::ostream * outputStream_default_ = &std::cout;
379 
380  //
381  // Current solver parameters and other values.
382  //
383 
385  MagnitudeType convtol_;
386 
388  MagnitudeType orthoKappa_;
389 
395  MagnitudeType achievedTol_;
396 
398  int maxIters_;
399 
401  int numIters_;
402 
404  int blockSize_, verbosity_, outputStyle_, outputFreq_;
405  bool adaptiveBlockSize_, showMaxResNormOnly_, useSingleReduction_;
406  std::string orthoType_, resScale_;
407  bool foldConvergenceDetectionIntoAllreduce_;
408 
410  std::string label_;
411 
413  Teuchos::RCP<Teuchos::Time> timerSolve_;
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  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
439  label_(label_default_),
440  isSet_(false)
441 {}
442 
443 
444 // Basic Constructor
445 template<class ScalarType, class MV, class OP>
449  problem_(problem),
450  outputStream_(Teuchos::rcp(outputStream_default_,false)),
451  convtol_(DefaultSolverParameters::convTol),
452  orthoKappa_(DefaultSolverParameters::orthoKappa),
453  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
454  maxIters_(maxIters_default_),
455  numIters_(0),
456  blockSize_(blockSize_default_),
457  verbosity_(verbosity_default_),
458  outputStyle_(outputStyle_default_),
459  outputFreq_(outputFreq_default_),
460  adaptiveBlockSize_(adaptiveBlockSize_default_),
461  showMaxResNormOnly_(showMaxResNormOnly_default_),
462  useSingleReduction_(useSingleReduction_default_),
463  orthoType_(orthoType_default_),
464  resScale_(resScale_default_),
465  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
466  label_(label_default_),
467  isSet_(false)
468 {
469  TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
470  "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
471 
472  // If the user passed in a nonnull parameter list, set parameters.
473  // Otherwise, the next solve() call will use default parameters,
474  // unless the user calls setParameters() first.
475  if (! pl.is_null()) {
476  setParameters (pl);
477  }
478 }
479 
480 template<class ScalarType, class MV, class OP>
481 void
484 {
485  // Create the internal parameter list if one doesn't already exist.
486  if (params_ == Teuchos::null) {
487  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
488  }
489  else {
490  params->validateParameters(*getValidParameters());
491  }
492 
493  // Check for maximum number of iterations
494  if (params->isParameter("Maximum Iterations")) {
495  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
496 
497  // Update parameter in our list and in status test.
498  params_->set("Maximum Iterations", maxIters_);
499  if (maxIterTest_!=Teuchos::null)
500  maxIterTest_->setMaxIters( maxIters_ );
501  }
502 
503  // Check for blocksize
504  if (params->isParameter("Block Size")) {
505  blockSize_ = params->get("Block Size",blockSize_default_);
506  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
507  "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
508 
509  // Update parameter in our list.
510  params_->set("Block Size", blockSize_);
511  }
512 
513  // Check if the blocksize should be adaptive
514  if (params->isParameter("Adaptive Block Size")) {
515  adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
516 
517  // Update parameter in our list.
518  params_->set("Adaptive Block Size", adaptiveBlockSize_);
519  }
520 
521  // Check if the user is requesting the single-reduction version of CG (only for blocksize == 1)
522  if (params->isParameter("Use Single Reduction")) {
523  useSingleReduction_ = params->get("Use Single Reduction", useSingleReduction_default_);
524  if (useSingleReduction_)
525  foldConvergenceDetectionIntoAllreduce_ = params->get("Fold Convergence Detection Into Allreduce",
526  foldConvergenceDetectionIntoAllreduce_default_);
527  }
528 
529  // Check to see if the timer label changed.
530  if (params->isParameter("Timer Label")) {
531  std::string tempLabel = params->get("Timer Label", label_default_);
532 
533  // Update parameter in our list and solver timer
534  if (tempLabel != label_) {
535  label_ = tempLabel;
536  params_->set("Timer Label", label_);
537  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
538 #ifdef BELOS_TEUCHOS_TIME_MONITOR
539  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
540 #endif
541  if (ortho_ != Teuchos::null) {
542  ortho_->setLabel( label_ );
543  }
544  }
545  }
546 
547  // Check if the orthogonalization changed.
548  if (params->isParameter("Orthogonalization")) {
549  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
550  TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS",
551  std::invalid_argument,
552  "Belos::BlockCGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
553  if (tempOrthoType != orthoType_) {
554  orthoType_ = tempOrthoType;
555  params_->set("Orthogonalization", orthoType_);
556  // Create orthogonalization manager
557  if (orthoType_=="DGKS") {
558  if (orthoKappa_ <= 0) {
559  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
560  }
561  else {
562  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
563  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
564  }
565  }
566  else if (orthoType_=="ICGS") {
567  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
568  }
569  else if (orthoType_=="IMGS") {
570  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
571  }
572  }
573  }
574 
575  // Check which orthogonalization constant to use.
576  if (params->isParameter("Orthogonalization Constant")) {
577  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
578  orthoKappa_ = params->get ("Orthogonalization Constant",
579  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
580  }
581  else {
582  orthoKappa_ = params->get ("Orthogonalization Constant",
584  }
585 
586  // Update parameter in our list.
587  params_->set("Orthogonalization Constant",orthoKappa_);
588  if (orthoType_=="DGKS") {
589  if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
590  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
591  }
592  }
593  }
594 
595  // Check for a change in verbosity level
596  if (params->isParameter("Verbosity")) {
597  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
598  verbosity_ = params->get("Verbosity", verbosity_default_);
599  } else {
600  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
601  }
602 
603  // Update parameter in our list.
604  params_->set("Verbosity", verbosity_);
605  if (printer_ != Teuchos::null)
606  printer_->setVerbosity(verbosity_);
607  }
608 
609  // Check for a change in output style
610  if (params->isParameter("Output Style")) {
611  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
612  outputStyle_ = params->get("Output Style", outputStyle_default_);
613  } else {
614  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
615  }
616 
617  // Update parameter in our list.
618  params_->set("Output Style", outputStyle_);
619  outputTest_ = Teuchos::null;
620  }
621 
622  // output stream
623  if (params->isParameter("Output Stream")) {
624  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
625 
626  // Update parameter in our list.
627  params_->set("Output Stream", outputStream_);
628  if (printer_ != Teuchos::null)
629  printer_->setOStream( outputStream_ );
630  }
631 
632  // frequency level
633  if (verbosity_ & Belos::StatusTestDetails) {
634  if (params->isParameter("Output Frequency")) {
635  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
636  }
637 
638  // Update parameter in out list and output status test.
639  params_->set("Output Frequency", outputFreq_);
640  if (outputTest_ != Teuchos::null)
641  outputTest_->setOutputFrequency( outputFreq_ );
642  }
643 
644  // Create output manager if we need to.
645  if (printer_ == Teuchos::null) {
646  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
647  }
648 
649  // Convergence
650  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
651  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
652 
653  // Check for convergence tolerance
654  if (params->isParameter("Convergence Tolerance")) {
655  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
656  convtol_ = params->get ("Convergence Tolerance",
657  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
658  }
659  else {
660  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
661  }
662 
663  // Update parameter in our list and residual tests.
664  params_->set("Convergence Tolerance", convtol_);
665  if (convTest_ != Teuchos::null)
666  convTest_->setTolerance( convtol_ );
667  }
668 
669  if (params->isParameter("Show Maximum Residual Norm Only")) {
670  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
671 
672  // Update parameter in our list and residual tests
673  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
674  if (convTest_ != Teuchos::null)
675  convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
676  }
677 
678  // Check for a change in scaling, if so we need to build new residual tests.
679  bool newResTest = false;
680  {
681  std::string tempResScale = resScale_;
682  if (params->isParameter ("Implicit Residual Scaling")) {
683  tempResScale = params->get<std::string> ("Implicit Residual Scaling");
684  }
685 
686  // Only update the scaling if it's different.
687  if (resScale_ != tempResScale) {
688  const Belos::ScaleType resScaleType =
689  convertStringToScaleType (tempResScale);
690  resScale_ = tempResScale;
691 
692  // Update parameter in our list and residual tests
693  params_->set ("Implicit Residual Scaling", resScale_);
694 
695  if (! convTest_.is_null ()) {
696  try {
697  NormType normType = Belos::TwoNorm;
698  if (params->isParameter("Residual Norm")) {
699  if (params->isType<std::string> ("Residual Norm")) {
700  normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
701  }
702  }
703  convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
704  convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
705  }
706  catch (std::exception& e) {
707  // Make sure the convergence test gets constructed again.
708  newResTest = true;
709  }
710  }
711  }
712  }
713 
714  // Create status tests if we need to.
715 
716  // Basic test checks maximum iterations and native residual.
717  if (maxIterTest_ == Teuchos::null)
718  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
719 
720  // Implicit residual test, using the native residual to determine if convergence was achieved.
721  if (convTest_.is_null () || newResTest) {
722 
723  NormType normType = Belos::TwoNorm;
724  if (params->isParameter("Residual Norm")) {
725  if (params->isType<std::string> ("Residual Norm")) {
726  normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
727  }
728  }
729 
730  convTest_ = rcp (new StatusTestResNorm_t (convtol_, 1, showMaxResNormOnly_));
731  convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
732  convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
733  }
734 
735  if (sTest_.is_null () || newResTest)
736  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
737 
738  if (outputTest_.is_null () || newResTest) {
739 
740  // Create the status test output class.
741  // This class manages and formats the output from the status test.
742  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
743  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
744 
745  // Set the solver string for the output test
746  std::string solverDesc = " Block CG ";
747  outputTest_->setSolverDesc( solverDesc );
748 
749  }
750 
751  // Create orthogonalization manager if we need to.
752  if (ortho_ == Teuchos::null) {
753  params_->set("Orthogonalization", orthoType_);
754  if (orthoType_=="DGKS") {
755  if (orthoKappa_ <= 0) {
756  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
757  }
758  else {
759  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
760  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
761  }
762  }
763  else if (orthoType_=="ICGS") {
764  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
765  }
766  else if (orthoType_=="IMGS") {
767  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
768  }
769  else {
770  TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
771  "Belos::BlockCGSolMgr(): Invalid orthogonalization type.");
772  }
773  }
774 
775  // Create the timer if we need to.
776  if (timerSolve_ == Teuchos::null) {
777  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
778 #ifdef BELOS_TEUCHOS_TIME_MONITOR
779  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
780 #endif
781  }
782 
783  // Inform the solver manager that the current parameters were set.
784  isSet_ = true;
785 }
786 
787 
788 template<class ScalarType, class MV, class OP>
791 {
793 
794  // Set all the valid parameters and their default values.
795  if(is_null(validPL)) {
796  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
797  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
798  "The relative residual tolerance that needs to be achieved by the\n"
799  "iterative solver in order for the linear system to be declared converged.");
800  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
801  "The maximum number of block iterations allowed for each\n"
802  "set of RHS solved.");
803  pl->set("Block Size", static_cast<int>(blockSize_default_),
804  "The number of vectors in each block.");
805  pl->set("Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
806  "Whether the solver manager should adapt to the block size\n"
807  "based on the number of RHS to solve.");
808  pl->set("Verbosity", static_cast<int>(verbosity_default_),
809  "What type(s) of solver information should be outputted\n"
810  "to the output stream.");
811  pl->set("Output Style", static_cast<int>(outputStyle_default_),
812  "What style is used for the solver information outputted\n"
813  "to the output stream.");
814  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
815  "How often convergence information should be outputted\n"
816  "to the output stream.");
817  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
818  "A reference-counted pointer to the output stream where all\n"
819  "solver output is sent.");
820  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
821  "When convergence information is printed, only show the maximum\n"
822  "relative residual norm when the block size is greater than one.");
823  pl->set("Use Single Reduction", static_cast<bool>(useSingleReduction_default_),
824  "Use single reduction iteration when the block size is one.");
825  pl->set("Implicit Residual Scaling", resScale_default_,
826  "The type of scaling used in the residual convergence test.");
827  pl->set("Timer Label", static_cast<const char *>(label_default_),
828  "The string to use as a prefix for the timer labels.");
829  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
830  "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
831  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
832  "The constant used by DGKS orthogonalization to determine\n"
833  "whether another step of classical Gram-Schmidt is necessary.");
834  pl->set("Residual Norm",static_cast<const char *>(resNorm_default_),
835  "Norm used for the convergence check on the residual.");
836  pl->set("Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
837  "Merge the allreduce for convergence detection with the one for CG.\n"
838  "This saves one all-reduce, but incurs more computation.");
839  validPL = pl;
840  }
841  return validPL;
842 }
843 
844 
845 // solve()
846 template<class ScalarType, class MV, class OP>
848  using Teuchos::RCP;
849  using Teuchos::rcp;
850  using Teuchos::rcp_const_cast;
851  using Teuchos::rcp_dynamic_cast;
852 
853  // Set the current parameters if they were not set before. NOTE:
854  // This may occur if the user generated the solver manager with the
855  // default constructor and then didn't set any parameters using
856  // setParameters().
857  if (!isSet_) {
858  setParameters(Teuchos::parameterList(*getValidParameters()));
859  }
860 
863 
864  TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(),
866  "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
867  "has not been called.");
868 
869  // Create indices for the linear systems to be solved.
870  int startPtr = 0;
871  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
872  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
873 
874  std::vector<int> currIdx, currIdx2;
875  // If an adaptive block size is allowed then only the linear
876  // systems that need to be solved are solved. Otherwise, the index
877  // set is generated that informs the linear problem that some
878  // linear systems are augmented.
879  if ( adaptiveBlockSize_ ) {
880  blockSize_ = numCurrRHS;
881  currIdx.resize( numCurrRHS );
882  currIdx2.resize( numCurrRHS );
883  for (int i=0; i<numCurrRHS; ++i)
884  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
885 
886  }
887  else {
888  currIdx.resize( blockSize_ );
889  currIdx2.resize( blockSize_ );
890  for (int i=0; i<numCurrRHS; ++i)
891  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
892  for (int i=numCurrRHS; i<blockSize_; ++i)
893  { currIdx[i] = -1; currIdx2[i] = i; }
894  }
895 
896  // Inform the linear problem of the current linear system to solve.
897  problem_->setLSIndex( currIdx );
898 
900  // Set up the parameter list for the Iteration subclass.
902  plist.set("Block Size",blockSize_);
903 
904  // Reset the output status test (controls all the other status tests).
905  outputTest_->reset();
906 
907  // Assume convergence is achieved, then let any failed convergence
908  // set this to false. "Innocent until proven guilty."
909  bool isConverged = true;
910 
912  // Set up the BlockCG Iteration subclass.
913 
914  RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
915  if (blockSize_ == 1) {
916  // Standard (nonblock) CG is faster for the special case of a
917  // block size of 1. A single reduction iteration can also be used
918  // if collectives are more expensive than vector updates.
919  if (useSingleReduction_) {
920  plist.set("Fold Convergence Detection Into Allreduce",
921  foldConvergenceDetectionIntoAllreduce_);
922  block_cg_iter =
923  rcp (new CGSingleRedIter<ScalarType,MV,OP> (problem_, printer_,
924  outputTest_, convTest_, plist));
925  }
926  else {
927  block_cg_iter =
928  rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_,
929  outputTest_, plist));
930  }
931  } else {
932  block_cg_iter =
933  rcp (new BlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_,
934  ortho_, plist));
935  }
936 
937 
938  // Enter solve() iterations
939  {
940 #ifdef BELOS_TEUCHOS_TIME_MONITOR
941  Teuchos::TimeMonitor slvtimer(*timerSolve_);
942 #endif
943 
944  while ( numRHS2Solve > 0 ) {
945  //
946  // Reset the active / converged vectors from this block
947  std::vector<int> convRHSIdx;
948  std::vector<int> currRHSIdx( currIdx );
949  currRHSIdx.resize(numCurrRHS);
950 
951  // Reset the number of iterations.
952  block_cg_iter->resetNumIters();
953 
954  // Reset the number of calls that the status test output knows about.
955  outputTest_->resetNumCalls();
956 
957  // Get the current residual for this block of linear systems.
958  RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
959 
960  // Set the new state and initialize the solver.
962  newstate.R = R_0;
963  block_cg_iter->initializeCG(newstate);
964 
965  while(1) {
966 
967  // tell block_cg_iter to iterate
968  try {
969  block_cg_iter->iterate();
970  //
971  // Check whether any of the linear systems converged.
972  //
973  if (convTest_->getStatus() == Passed) {
974  // At least one of the linear system(s) converged.
975  //
976  // Get the column indices of the linear systems that converged.
977  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
978  std::vector<int> convIdx =
979  rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
980 
981  // If the number of converged linear systems equals the
982  // number of linear systems currently being solved, then
983  // we are done with this block.
984  if (convIdx.size() == currRHSIdx.size())
985  break; // break from while(1){block_cg_iter->iterate()}
986 
987  // Inform the linear problem that we are finished with
988  // this current linear system.
989  problem_->setCurrLS();
990 
991  // Reset currRHSIdx to contain the right-hand sides that
992  // are left to converge for this block.
993  int have = 0;
994  std::vector<int> unconvIdx(currRHSIdx.size());
995  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
996  bool found = false;
997  for (unsigned int j=0; j<convIdx.size(); ++j) {
998  if (currRHSIdx[i] == convIdx[j]) {
999  found = true;
1000  break;
1001  }
1002  }
1003  if (!found) {
1004  currIdx2[have] = currIdx2[i];
1005  currRHSIdx[have++] = currRHSIdx[i];
1006  }
1007  else {
1008  }
1009  }
1010  currRHSIdx.resize(have);
1011  currIdx2.resize(have);
1012 
1013  // Set the remaining indices after deflation.
1014  problem_->setLSIndex( currRHSIdx );
1015 
1016  // Get the current residual vector.
1017  std::vector<MagnitudeType> norms;
1018  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
1019  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
1020 
1021  // Set the new blocksize for the solver.
1022  block_cg_iter->setBlockSize( have );
1023 
1024  // Set the new state and initialize the solver.
1026  defstate.R = R_0;
1027  block_cg_iter->initializeCG(defstate);
1028  }
1029  //
1030  // None of the linear systems converged. Check whether the
1031  // maximum iteration count was reached.
1032  //
1033  else if (maxIterTest_->getStatus() == Passed) {
1034  isConverged = false; // None of the linear systems converged.
1035  break; // break from while(1){block_cg_iter->iterate()}
1036  }
1037  //
1038  // iterate() returned, but none of our status tests Passed.
1039  // This indicates a bug.
1040  //
1041  else {
1042  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1043  "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
1044  "the maximum iteration count test passed. Please report this bug "
1045  "to the Belos developers.");
1046  }
1047  }
1048  catch (const std::exception &e) {
1049  std::ostream& err = printer_->stream (Errors);
1050  err << "Error! Caught std::exception in CGIteration::iterate() at "
1051  << "iteration " << block_cg_iter->getNumIters() << std::endl
1052  << e.what() << std::endl;
1053  throw;
1054  }
1055  }
1056 
1057  // Inform the linear problem that we are finished with this
1058  // block linear system.
1059  problem_->setCurrLS();
1060 
1061  // Update indices for the linear systems to be solved.
1062  startPtr += numCurrRHS;
1063  numRHS2Solve -= numCurrRHS;
1064  if ( numRHS2Solve > 0 ) {
1065  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1066 
1067  if ( adaptiveBlockSize_ ) {
1068  blockSize_ = numCurrRHS;
1069  currIdx.resize( numCurrRHS );
1070  currIdx2.resize( numCurrRHS );
1071  for (int i=0; i<numCurrRHS; ++i)
1072  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1073  }
1074  else {
1075  currIdx.resize( blockSize_ );
1076  currIdx2.resize( blockSize_ );
1077  for (int i=0; i<numCurrRHS; ++i)
1078  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1079  for (int i=numCurrRHS; i<blockSize_; ++i)
1080  { currIdx[i] = -1; currIdx2[i] = i; }
1081  }
1082  // Set the next indices.
1083  problem_->setLSIndex( currIdx );
1084 
1085  // Set the new blocksize for the solver.
1086  block_cg_iter->setBlockSize( blockSize_ );
1087  }
1088  else {
1089  currIdx.resize( numRHS2Solve );
1090  }
1091 
1092  }// while ( numRHS2Solve > 0 )
1093 
1094  }
1095 
1096  // print final summary
1097  sTest_->print( printer_->stream(FinalSummary) );
1098 
1099  // print timing information
1100 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1101  // Calling summarize() requires communication in general, so don't
1102  // call it unless the user wants to print out timing details.
1103  // summarize() will do all the work even if it's passed a "black
1104  // hole" output stream.
1105  if (verbosity_ & TimingDetails) {
1106  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1107  }
1108 #endif
1109 
1110  // Save the iteration count for this solve.
1111  numIters_ = maxIterTest_->getNumIters();
1112 
1113  // Save the convergence test value ("achieved tolerance") for this solve.
1114  {
1115  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1116  // testValues is nonnull and not persistent.
1117  const std::vector<MagnitudeType>* pTestValues =
1118  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1119 
1120  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1121  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1122  "method returned NULL. Please report this bug to the Belos developers.");
1123 
1124  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1125  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1126  "method returned a vector of length zero. Please report this bug to the "
1127  "Belos developers.");
1128 
1129  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1130  // achieved tolerances for all vectors in the current solve(), or
1131  // just for the vectors from the last deflation?
1132  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1133  }
1134 
1135  if (!isConverged) {
1136  return Unconverged; // return from BlockCGSolMgr::solve()
1137  }
1138  return Converged; // return from BlockCGSolMgr::solve()
1139 }
1140 
1141 // This method requires the solver manager to return a std::string that describes itself.
1142 template<class ScalarType, class MV, class OP>
1144 {
1145  std::ostringstream oss;
1146  oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1147  oss << "{";
1148  oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
1149  oss << "}";
1150  return oss.str();
1151 }
1152 
1153 } // end Belos namespace
1154 
1155 #endif /* BELOS_BLOCK_CG_SOLMGR_HPP */
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:307
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...
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp: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
Structure to contain pointers to CGIteration state variables.
T & get(const std::string &name, T def_value)
Belos concrete class for performing the conjugate-gradient (CG) iteration.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
Belos concrete class for performing the block conjugate-gradient (CG) iteration.
An implementation of StatusTestResNorm using a family of residual norms.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
The Belos::BlockCGSolMgr provides a powerful and fully-featured solver manager over the CG and BlockC...
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:301
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.
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)
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
bool isParameter(const std::string &name) const
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: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)
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.
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...
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
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 ...
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.
bool isType(const std::string &name) const
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:291
This class implements the preconditioned single-reduction Conjugate Gradient (CG) iteration...
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...
bool is_null() const

Generated on Fri Jun 5 2020 10:20:37 for Belos by doxygen 1.8.5