Belos  Version of the Day
BelosBlockCGSolMgr.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_BLOCK_CG_SOLMGR_HPP
43 #define BELOS_BLOCK_CG_SOLMGR_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
55 #include "BelosCGIter.hpp"
56 #include "BelosCGSingleRedIter.hpp"
57 #include "BelosBlockCGIter.hpp"
61 #include "BelosStatusTestCombo.hpp"
63 #include "BelosOutputManager.hpp"
64 #include "Teuchos_LAPACK.hpp"
65 #ifdef BELOS_TEUCHOS_TIME_MONITOR
66 # include "Teuchos_TimeMonitor.hpp"
67 #endif
68 #if defined(HAVE_TEUCHOSCORE_CXX11)
69 # include <type_traits>
70 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
71 #include <algorithm>
72 
89 namespace Belos {
90 
92 
93 
101  BlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
102  {}};
103 
104  template<class ScalarType, class MV, class OP,
105  const bool lapackSupportsScalarType =
108  public Details::SolverManagerRequiresLapack<ScalarType,MV,OP>
109  {
110  static const bool requiresLapack =
112  typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
113  requiresLapack> base_type;
114 
115  public:
117  base_type ()
118  {}
119  BlockCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem,
120  const Teuchos::RCP<Teuchos::ParameterList>& pl) :
121  base_type ()
122  {}
123  virtual ~BlockCGSolMgr () {}
124  };
125 
126 
127  // Partial specialization for ScalarType types for which
128  // Teuchos::LAPACK has a valid implementation. This contains the
129  // actual working implementation of BlockCGSolMgr.
130  template<class ScalarType, class MV, class OP>
131  class BlockCGSolMgr<ScalarType, MV, OP, true> :
132  public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
133  {
134  // This partial specialization is already chosen for those scalar types
135  // that support lapack, so we don't need to have an additional compile-time
136  // check that the scalar type is float/double/complex.
137 // #if defined(HAVE_TEUCHOSCORE_CXX11)
138 // # if defined(HAVE_TEUCHOS_COMPLEX)
139 // static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
140 // std::is_same<ScalarType, std::complex<double> >::value ||
141 // std::is_same<ScalarType, float>::value ||
142 // std::is_same<ScalarType, double>::value,
143 // "Belos::GCRODRSolMgr: ScalarType must be one of the four "
144 // "types (S,D,C,Z) supported by LAPACK.");
145 // # else
146 // static_assert (std::is_same<ScalarType, float>::value ||
147 // std::is_same<ScalarType, double>::value,
148 // "Belos::GCRODRSolMgr: ScalarType must be float or double. "
149 // "Complex arithmetic support is currently disabled. To "
150 // "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
151 // # endif // defined(HAVE_TEUCHOS_COMPLEX)
152 // #endif // defined(HAVE_TEUCHOSCORE_CXX11)
153 
154  private:
157  typedef Teuchos::ScalarTraits<ScalarType> SCT;
158  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
159  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
160 
161  public:
162 
164 
165 
171  BlockCGSolMgr();
172 
209  BlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
210  const Teuchos::RCP<Teuchos::ParameterList> &pl );
211 
213  virtual ~BlockCGSolMgr() {};
214 
216  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
217  return Teuchos::rcp(new BlockCGSolMgr<ScalarType,MV,OP>);
218  }
220 
222 
223 
224  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
225  return *problem_;
226  }
227 
230  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
231 
234  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
235 
241  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
242  return Teuchos::tuple(timerSolve_);
243  }
244 
250  MagnitudeType achievedTol() const override {
251  return achievedTol_;
252  }
253 
255  int getNumIters() const override {
256  return numIters_;
257  }
258 
261  bool isLOADetected() const override { return false; }
262 
264 
266 
267 
269  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
270 
272  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
273 
275 
277 
278 
282  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
284 
286 
287 
305  ReturnType solve() override;
306 
308 
311 
313  std::string description() const override;
314 
316 
317  private:
318 
320  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
321 
323  Teuchos::RCP<OutputManager<ScalarType> > printer_;
325  Teuchos::RCP<std::ostream> outputStream_;
326 
331  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
332 
334  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
335 
337  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
338 
340  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
341 
343  Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
344 
346  Teuchos::RCP<Teuchos::ParameterList> params_;
347 
348  //
349  // Default solver parameters.
350  //
351  static constexpr int maxIters_default_ = 1000;
352  static constexpr bool adaptiveBlockSize_default_ = true;
353  static constexpr bool showMaxResNormOnly_default_ = false;
354  static constexpr bool useSingleReduction_default_ = false;
355  static constexpr int blockSize_default_ = 1;
356  static constexpr int verbosity_default_ = Belos::Errors;
357  static constexpr int outputStyle_default_ = Belos::General;
358  static constexpr int outputFreq_default_ = -1;
359  static constexpr const char * resNorm_default_ = "TwoNorm";
360  static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ = false;
361  static constexpr const char * resScale_default_ = "Norm of Initial Residual";
362  static constexpr const char * label_default_ = "Belos";
363  static constexpr const char * orthoType_default_ = "ICGS";
364  static constexpr bool assertPositiveDefiniteness_default_ = true;
365 // https://stackoverflow.com/questions/24398102/constexpr-and-initialization-of-a-static-const-void-pointer-with-reinterpret-cas
366 #if defined(_WIN32) && defined(__clang__)
367  static constexpr std::ostream * outputStream_default_ =
368  __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
369 #else
370  static constexpr std::ostream * outputStream_default_ = &std::cout;
371 #endif
372 
373  //
374  // Current solver parameters and other values.
375  //
376 
378  MagnitudeType convtol_;
379 
381  MagnitudeType orthoKappa_;
382 
388  MagnitudeType achievedTol_;
389 
391  int maxIters_;
392 
394  int numIters_;
395 
397  int blockSize_, verbosity_, outputStyle_, outputFreq_;
398  bool adaptiveBlockSize_, showMaxResNormOnly_, useSingleReduction_;
399  std::string orthoType_, resScale_;
400  bool assertPositiveDefiniteness_;
401  bool foldConvergenceDetectionIntoAllreduce_;
402 
404  std::string label_;
405 
407  Teuchos::RCP<Teuchos::Time> timerSolve_;
408 
410  bool isSet_;
411  };
412 
413 
414 // Empty Constructor
415 template<class ScalarType, class MV, class OP>
417  outputStream_(Teuchos::rcp(outputStream_default_,false)),
418  convtol_(DefaultSolverParameters::convTol),
419  orthoKappa_(DefaultSolverParameters::orthoKappa),
420  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
421  maxIters_(maxIters_default_),
422  numIters_(0),
423  blockSize_(blockSize_default_),
424  verbosity_(verbosity_default_),
425  outputStyle_(outputStyle_default_),
426  outputFreq_(outputFreq_default_),
427  adaptiveBlockSize_(adaptiveBlockSize_default_),
428  showMaxResNormOnly_(showMaxResNormOnly_default_),
429  useSingleReduction_(useSingleReduction_default_),
430  orthoType_(orthoType_default_),
431  resScale_(resScale_default_),
432  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
433  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
434  label_(label_default_),
435  isSet_(false)
436 {}
437 
438 
439 // Basic Constructor
440 template<class ScalarType, class MV, class OP>
442 BlockCGSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
443  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
444  problem_(problem),
445  outputStream_(Teuchos::rcp(outputStream_default_,false)),
446  convtol_(DefaultSolverParameters::convTol),
447  orthoKappa_(DefaultSolverParameters::orthoKappa),
448  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
449  maxIters_(maxIters_default_),
450  numIters_(0),
451  blockSize_(blockSize_default_),
452  verbosity_(verbosity_default_),
453  outputStyle_(outputStyle_default_),
454  outputFreq_(outputFreq_default_),
455  adaptiveBlockSize_(adaptiveBlockSize_default_),
456  showMaxResNormOnly_(showMaxResNormOnly_default_),
457  useSingleReduction_(useSingleReduction_default_),
458  orthoType_(orthoType_default_),
459  resScale_(resScale_default_),
460  assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
461  foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
462  label_(label_default_),
463  isSet_(false)
464 {
465  TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
466  "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
467 
468  // If the user passed in a nonnull parameter list, set parameters.
469  // Otherwise, the next solve() call will use default parameters,
470  // unless the user calls setParameters() first.
471  if (! pl.is_null()) {
472  setParameters (pl);
473  }
474 }
475 
476 template<class ScalarType, class MV, class OP>
477 void
479 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
480 {
481  // Create the internal parameter list if one doesn't already exist.
482  if (params_ == Teuchos::null) {
483  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
484  }
485  else {
486  params->validateParameters(*getValidParameters());
487  }
488 
489  // Check for maximum number of iterations
490  if (params->isParameter("Maximum Iterations")) {
491  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
492 
493  // Update parameter in our list and in status test.
494  params_->set("Maximum Iterations", maxIters_);
495  if (maxIterTest_!=Teuchos::null)
496  maxIterTest_->setMaxIters( maxIters_ );
497  }
498 
499  // Check for blocksize
500  if (params->isParameter("Block Size")) {
501  blockSize_ = params->get("Block Size",blockSize_default_);
502  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
503  "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
504 
505  // Update parameter in our list.
506  params_->set("Block Size", blockSize_);
507  }
508 
509  // Check if the blocksize should be adaptive
510  if (params->isParameter("Adaptive Block Size")) {
511  adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
512 
513  // Update parameter in our list.
514  params_->set("Adaptive Block Size", adaptiveBlockSize_);
515  }
516 
517  // Check if the user is requesting the single-reduction version of CG (only for blocksize == 1)
518  if (params->isParameter("Use Single Reduction")) {
519  useSingleReduction_ = params->get("Use Single Reduction", useSingleReduction_default_);
520  }
521 
522  if (params->isParameter("Fold Convergence Detection Into Allreduce")) {
523  foldConvergenceDetectionIntoAllreduce_ = params->get("Fold Convergence Detection Into Allreduce",
524  foldConvergenceDetectionIntoAllreduce_default_);
525  }
526 
527  // Check to see if the timer label changed.
528  if (params->isParameter("Timer Label")) {
529  std::string tempLabel = params->get("Timer Label", label_default_);
530 
531  // Update parameter in our list and solver timer
532  if (tempLabel != label_) {
533  label_ = tempLabel;
534  params_->set("Timer Label", label_);
535  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
536 #ifdef BELOS_TEUCHOS_TIME_MONITOR
537  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
538 #endif
539  if (ortho_ != Teuchos::null) {
540  ortho_->setLabel( label_ );
541  }
542  }
543  }
544 
545  // Check for a change in verbosity level
546  if (params->isParameter("Verbosity")) {
547  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
548  verbosity_ = params->get("Verbosity", verbosity_default_);
549  } else {
550  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
551  }
552 
553  // Update parameter in our list.
554  params_->set("Verbosity", verbosity_);
555  if (printer_ != Teuchos::null)
556  printer_->setVerbosity(verbosity_);
557  }
558 
559  // Check for a change in output style
560  if (params->isParameter("Output Style")) {
561  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
562  outputStyle_ = params->get("Output Style", outputStyle_default_);
563  } else {
564  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
565  }
566 
567  // Update parameter in our list.
568  params_->set("Output Style", outputStyle_);
569  outputTest_ = Teuchos::null;
570  }
571 
572  // output stream
573  if (params->isParameter("Output Stream")) {
574  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
575 
576  // Update parameter in our list.
577  params_->set("Output Stream", outputStream_);
578  if (printer_ != Teuchos::null)
579  printer_->setOStream( outputStream_ );
580  }
581 
582  // frequency level
583  if (verbosity_ & Belos::StatusTestDetails) {
584  if (params->isParameter("Output Frequency")) {
585  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
586  }
587 
588  // Update parameter in out list and output status test.
589  params_->set("Output Frequency", outputFreq_);
590  if (outputTest_ != Teuchos::null)
591  outputTest_->setOutputFrequency( outputFreq_ );
592  }
593 
594  // Create output manager if we need to.
595  if (printer_ == Teuchos::null) {
596  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
597  }
598 
599  // Check if the orthogonalization changed.
600  bool changedOrthoType = false;
601  if (params->isParameter("Orthogonalization")) {
602  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
603  if (tempOrthoType != orthoType_) {
604  orthoType_ = tempOrthoType;
605  changedOrthoType = true;
606  }
607  }
608  params_->set("Orthogonalization", orthoType_);
609 
610  // Check which orthogonalization constant to use.
611  if (params->isParameter("Orthogonalization Constant")) {
612  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
613  orthoKappa_ = params->get ("Orthogonalization Constant",
614  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
615  }
616  else {
617  orthoKappa_ = params->get ("Orthogonalization Constant",
619  }
620 
621  // Update parameter in our list.
622  params_->set("Orthogonalization Constant",orthoKappa_);
623  if (orthoType_=="DGKS") {
624  if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
625  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
626  }
627  }
628  }
629 
630  // Create orthogonalization manager if we need to.
631  if (ortho_ == Teuchos::null || changedOrthoType) {
633  Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
634  if (orthoType_=="DGKS" && orthoKappa_ > 0) {
635  paramsOrtho->set ("depTol", orthoKappa_ );
636  }
637 
638  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
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  NormType normType = Belos::TwoNorm;
690  if (params->isParameter("Residual Norm")) {
691  if (params->isType<std::string> ("Residual Norm")) {
692  normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
693  }
694  }
695  convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
696  convTest_->defineScaleForm (resScaleType, Belos::TwoNorm);
697  }
698  catch (std::exception& e) {
699  // Make sure the convergence test gets constructed again.
700  newResTest = true;
701  }
702  }
703  }
704  }
705 
706  // Create status tests if we need to.
707 
708  // Basic test checks maximum iterations and native residual.
709  if (maxIterTest_ == Teuchos::null)
710  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
711 
712  // Implicit residual test, using the native residual to determine if convergence was achieved.
713  if (convTest_.is_null () || newResTest) {
714 
715  NormType normType = Belos::TwoNorm;
716  if (params->isParameter("Residual Norm")) {
717  if (params->isType<std::string> ("Residual Norm")) {
718  normType = convertStringToNormType(params->get<std::string> ("Residual Norm"));
719  }
720  }
721 
722  convTest_ = rcp (new StatusTestResNorm_t (convtol_, 1, showMaxResNormOnly_));
723  convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
724  convTest_->defineScaleForm (convertStringToScaleType (resScale_), Belos::TwoNorm);
725  }
726 
727  if (sTest_.is_null () || newResTest)
728  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
729 
730  if (outputTest_.is_null () || newResTest) {
731 
732  // Create the status test output class.
733  // This class manages and formats the output from the status test.
734  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
735  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
736 
737  // Set the solver string for the output test
738  std::string solverDesc = " Block CG ";
739  outputTest_->setSolverDesc( solverDesc );
740 
741  }
742 
743  // BelosCgIter accepts a parameter specifying whether to assert for the positivity of p^H*A*p in the CG iteration
744  if (params->isParameter("Assert Positive Definiteness")) {
745  assertPositiveDefiniteness_ = Teuchos::getParameter<bool>(*params,"Assert Positive Definiteness");
746  params_->set("Assert Positive Definiteness", assertPositiveDefiniteness_);
747  }
748 
749  // Create the timer if we need to.
750  if (timerSolve_ == Teuchos::null) {
751  std::string solveLabel = label_ + ": BlockCGSolMgr total solve time";
752 #ifdef BELOS_TEUCHOS_TIME_MONITOR
753  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
754 #endif
755  }
756 
757  // Inform the solver manager that the current parameters were set.
758  isSet_ = true;
759 }
760 
761 
762 template<class ScalarType, class MV, class OP>
763 Teuchos::RCP<const Teuchos::ParameterList>
765 {
766  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
767 
768  // Set all the valid parameters and their default values.
769  if(is_null(validPL)) {
770  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
771  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
772  "The relative residual tolerance that needs to be achieved by the\n"
773  "iterative solver in order for the linear system to be declared converged.");
774  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
775  "The maximum number of block iterations allowed for each\n"
776  "set of RHS solved.");
777  pl->set("Block Size", static_cast<int>(blockSize_default_),
778  "The number of vectors in each block.");
779  pl->set("Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
780  "Whether the solver manager should adapt to the block size\n"
781  "based on the number of RHS to solve.");
782  pl->set("Verbosity", static_cast<int>(verbosity_default_),
783  "What type(s) of solver information should be outputted\n"
784  "to the output stream.");
785  pl->set("Output Style", static_cast<int>(outputStyle_default_),
786  "What style is used for the solver information outputted\n"
787  "to the output stream.");
788  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
789  "How often convergence information should be outputted\n"
790  "to the output stream.");
791  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
792  "A reference-counted pointer to the output stream where all\n"
793  "solver output is sent.");
794  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
795  "When convergence information is printed, only show the maximum\n"
796  "relative residual norm when the block size is greater than one.");
797  pl->set("Use Single Reduction", static_cast<bool>(useSingleReduction_default_),
798  "Use single reduction iteration when the block size is one.");
799  pl->set("Implicit Residual Scaling", resScale_default_,
800  "The type of scaling used in the residual convergence test.");
801  pl->set("Timer Label", static_cast<const char *>(label_default_),
802  "The string to use as a prefix for the timer labels.");
803  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
804  "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
805  pl->set("Assert Positive Definiteness",static_cast<bool>(assertPositiveDefiniteness_default_),
806  "Assert for positivity of p^H*A*p in CG iteration.");
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  pl->set("Residual Norm",static_cast<const char *>(resNorm_default_),
811  "Norm used for the convergence check on the residual.");
812  pl->set("Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
813  "Merge the allreduce for convergence detection with the one for CG.\n"
814  "This saves one all-reduce, but incurs more computation.");
815  validPL = pl;
816  }
817  return validPL;
818 }
819 
820 
821 // solve()
822 template<class ScalarType, class MV, class OP>
824  using Teuchos::RCP;
825  using Teuchos::rcp;
826  using Teuchos::rcp_const_cast;
827  using Teuchos::rcp_dynamic_cast;
828 
829  // Set the current parameters if they were not set before. NOTE:
830  // This may occur if the user generated the solver manager with the
831  // default constructor and then didn't set any parameters using
832  // setParameters().
833  if (!isSet_) {
834  setParameters(Teuchos::parameterList(*getValidParameters()));
835  }
836 
837  Teuchos::LAPACK<int,ScalarType> lapack;
838 
839  TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(),
841  "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
842  "has not been called.");
843 
844  // Create indices for the linear systems to be solved.
845  int startPtr = 0;
846  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
847  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
848 
849  std::vector<int> currIdx, currIdx2;
850  // If an adaptive block size is allowed then only the linear
851  // systems that need to be solved are solved. Otherwise, the index
852  // set is generated that informs the linear problem that some
853  // linear systems are augmented.
854  if ( adaptiveBlockSize_ ) {
855  blockSize_ = numCurrRHS;
856  currIdx.resize( numCurrRHS );
857  currIdx2.resize( numCurrRHS );
858  for (int i=0; i<numCurrRHS; ++i)
859  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
860 
861  }
862  else {
863  currIdx.resize( blockSize_ );
864  currIdx2.resize( blockSize_ );
865  for (int i=0; i<numCurrRHS; ++i)
866  { currIdx[i] = startPtr+i; currIdx2[i]=i; }
867  for (int i=numCurrRHS; i<blockSize_; ++i)
868  { currIdx[i] = -1; currIdx2[i] = i; }
869  }
870 
871  // Inform the linear problem of the current linear system to solve.
872  problem_->setLSIndex( currIdx );
873 
875  // Set up the parameter list for the Iteration subclass.
876  Teuchos::ParameterList plist;
877  plist.set("Block Size",blockSize_);
878 
879  // Reset the output status test (controls all the other status tests).
880  outputTest_->reset();
881 
882  // Assume convergence is achieved, then let any failed convergence
883  // set this to false. "Innocent until proven guilty."
884  bool isConverged = true;
885 
887  // Set up the BlockCG Iteration subclass.
888 
889  plist.set("Assert Positive Definiteness", assertPositiveDefiniteness_);
890 
891  RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
892  if (blockSize_ == 1) {
893  // Standard (nonblock) CG is faster for the special case of a
894  // block size of 1. A single reduction iteration can also be used
895  // if collectives are more expensive than vector updates.
896  plist.set("Fold Convergence Detection Into Allreduce",
897  foldConvergenceDetectionIntoAllreduce_);
898  if (useSingleReduction_) {
899  block_cg_iter =
900  rcp (new CGSingleRedIter<ScalarType,MV,OP> (problem_, printer_,
901  outputTest_, convTest_, plist));
902  }
903  else {
904  block_cg_iter =
905  rcp (new CGIter<ScalarType,MV,OP> (problem_, printer_,
906  outputTest_, convTest_, plist));
907  }
908  } else {
909  block_cg_iter =
910  rcp (new BlockCGIter<ScalarType,MV,OP> (problem_, printer_, outputTest_,
911  ortho_, plist));
912  }
913 
914 
915  // Enter solve() iterations
916  {
917 #ifdef BELOS_TEUCHOS_TIME_MONITOR
918  Teuchos::TimeMonitor slvtimer(*timerSolve_);
919 #endif
920 
921  while ( numRHS2Solve > 0 ) {
922  //
923  // Reset the active / converged vectors from this block
924  std::vector<int> convRHSIdx;
925  std::vector<int> currRHSIdx( currIdx );
926  currRHSIdx.resize(numCurrRHS);
927 
928  // Reset the number of iterations.
929  block_cg_iter->resetNumIters();
930 
931  // Reset the number of calls that the status test output knows about.
932  outputTest_->resetNumCalls();
933 
934  // Get the current residual for this block of linear systems.
935  RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
936 
937  // Set the new state and initialize the solver.
939  newstate.R = R_0;
940  block_cg_iter->initializeCG(newstate);
941 
942  while(1) {
943 
944  // tell block_cg_iter to iterate
945  try {
946  block_cg_iter->iterate();
947  //
948  // Check whether any of the linear systems converged.
949  //
950  if (convTest_->getStatus() == Passed) {
951  // At least one of the linear system(s) converged.
952  //
953  // Get the column indices of the linear systems that converged.
954  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
955  std::vector<int> convIdx =
956  rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
957 
958  // If the number of converged linear systems equals the
959  // number of linear systems currently being solved, then
960  // we are done with this block.
961  if (convIdx.size() == currRHSIdx.size())
962  break; // break from while(1){block_cg_iter->iterate()}
963 
964  // Inform the linear problem that we are finished with
965  // this current linear system.
966  problem_->setCurrLS();
967 
968  // Reset currRHSIdx to contain the right-hand sides that
969  // are left to converge for this block.
970  int have = 0;
971  std::vector<int> unconvIdx(currRHSIdx.size());
972  for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
973  bool found = false;
974  for (unsigned int j=0; j<convIdx.size(); ++j) {
975  if (currRHSIdx[i] == convIdx[j]) {
976  found = true;
977  break;
978  }
979  }
980  if (!found) {
981  currIdx2[have] = currIdx2[i];
982  currRHSIdx[have++] = currRHSIdx[i];
983  }
984  else {
985  }
986  }
987  currRHSIdx.resize(have);
988  currIdx2.resize(have);
989 
990  // Set the remaining indices after deflation.
991  problem_->setLSIndex( currRHSIdx );
992 
993  // Get the current residual vector.
994  std::vector<MagnitudeType> norms;
995  R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
996  for (int i=0; i<have; ++i) { currIdx2[i] = i; }
997 
998  // Set the new blocksize for the solver.
999  block_cg_iter->setBlockSize( have );
1000 
1001  // Set the new state and initialize the solver.
1003  defstate.R = R_0;
1004  block_cg_iter->initializeCG(defstate);
1005  }
1006  //
1007  // None of the linear systems converged. Check whether the
1008  // maximum iteration count was reached.
1009  //
1010  else if (maxIterTest_->getStatus() == Passed) {
1011  isConverged = false; // None of the linear systems converged.
1012  break; // break from while(1){block_cg_iter->iterate()}
1013  }
1014  //
1015  // iterate() returned, but none of our status tests Passed.
1016  // This indicates a bug.
1017  //
1018  else {
1019  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1020  "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
1021  "the maximum iteration count test passed. Please report this bug "
1022  "to the Belos developers.");
1023  }
1024  }
1025  catch (const std::exception &e) {
1026  std::ostream& err = printer_->stream (Errors);
1027  err << "Error! Caught std::exception in CGIteration::iterate() at "
1028  << "iteration " << block_cg_iter->getNumIters() << std::endl
1029  << e.what() << std::endl;
1030  throw;
1031  }
1032  }
1033 
1034  // Inform the linear problem that we are finished with this
1035  // block linear system.
1036  problem_->setCurrLS();
1037 
1038  // Update indices for the linear systems to be solved.
1039  startPtr += numCurrRHS;
1040  numRHS2Solve -= numCurrRHS;
1041  if ( numRHS2Solve > 0 ) {
1042  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1043 
1044  if ( adaptiveBlockSize_ ) {
1045  blockSize_ = numCurrRHS;
1046  currIdx.resize( numCurrRHS );
1047  currIdx2.resize( numCurrRHS );
1048  for (int i=0; i<numCurrRHS; ++i)
1049  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1050  }
1051  else {
1052  currIdx.resize( blockSize_ );
1053  currIdx2.resize( blockSize_ );
1054  for (int i=0; i<numCurrRHS; ++i)
1055  { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1056  for (int i=numCurrRHS; i<blockSize_; ++i)
1057  { currIdx[i] = -1; currIdx2[i] = i; }
1058  }
1059  // Set the next indices.
1060  problem_->setLSIndex( currIdx );
1061 
1062  // Set the new blocksize for the solver.
1063  block_cg_iter->setBlockSize( blockSize_ );
1064  }
1065  else {
1066  currIdx.resize( numRHS2Solve );
1067  }
1068 
1069  }// while ( numRHS2Solve > 0 )
1070 
1071  }
1072 
1073  // print final summary
1074  sTest_->print( printer_->stream(FinalSummary) );
1075 
1076  // print timing information
1077 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1078  // Calling summarize() requires communication in general, so don't
1079  // call it unless the user wants to print out timing details.
1080  // summarize() will do all the work even if it's passed a "black
1081  // hole" output stream.
1082  if (verbosity_ & TimingDetails) {
1083  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1084  }
1085 #endif
1086 
1087  // Save the iteration count for this solve.
1088  numIters_ = maxIterTest_->getNumIters();
1089 
1090  // Save the convergence test value ("achieved tolerance") for this solve.
1091  {
1092  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1093  // testValues is nonnull and not persistent.
1094  const std::vector<MagnitudeType>* pTestValues =
1095  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1096 
1097  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1098  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1099  "method returned NULL. Please report this bug to the Belos developers.");
1100 
1101  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1102  "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1103  "method returned a vector of length zero. Please report this bug to the "
1104  "Belos developers.");
1105 
1106  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1107  // achieved tolerances for all vectors in the current solve(), or
1108  // just for the vectors from the last deflation?
1109  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1110  }
1111 
1112  if (!isConverged) {
1113  return Unconverged; // return from BlockCGSolMgr::solve()
1114  }
1115  return Converged; // return from BlockCGSolMgr::solve()
1116 }
1117 
1118 // This method requires the solver manager to return a std::string that describes itself.
1119 template<class ScalarType, class MV, class OP>
1121 {
1122  std::ostringstream oss;
1123  oss << "Belos::BlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1124  oss << "{";
1125  oss << "Ortho Type='"<<orthoType_<<"\', Block Size=" << blockSize_;
1126  oss << "}";
1127  return oss.str();
1128 }
1129 
1130 } // end Belos namespace
1131 
1132 #endif /* BELOS_BLOCK_CG_SOLMGR_HPP */
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:299
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
Teuchos::RCP< const MV > R
The current residual.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Collection of types and exceptions used within the Belos solvers.
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.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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:79
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.
The Belos::BlockCGSolMgr provides a powerful and fully-featured solver manager over the CG and BlockC...
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:293
Belos::StatusTest class for specifying a maximum number of iterations.
BlockCGSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
BlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
A factory class for generating StatusTestOutput objects.
BlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
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...
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.
NormType convertStringToNormType(const std::string &normType)
Convert the given string to its NormType enum value.
Definition: BelosTypes.cpp:88
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
MagnitudeType achievedTol() 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().
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
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.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:283
This class implements the preconditioned single-reduction Conjugate Gradient (CG) iteration...
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...
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
Return an instance of the specified MatOrthoManager subclass.

Generated for Belos by doxygen 1.8.14