Belos Package Browser (Single Doxygen Collection)  Development
BelosBlockGmresSolMgr.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_GMRES_SOLMGR_HPP
43 #define BELOS_BLOCK_GMRES_SOLMGR_HPP
44 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
55 #include "BelosGmresIteration.hpp"
56 #include "BelosBlockGmresIter.hpp"
57 #include "BelosBlockFGmresIter.hpp"
62 #include "BelosStatusTestCombo.hpp"
65 #include "BelosOutputManager.hpp"
66 #ifdef BELOS_TEUCHOS_TIME_MONITOR
67 #include "Teuchos_TimeMonitor.hpp"
68 #endif
69 
80 namespace Belos {
81 
83 
84 
92  BlockGmresSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
93  {}};
94 
102  BlockGmresSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
103  {}};
104 
121 template<class ScalarType, class MV, class OP>
122 class BlockGmresSolMgr : public SolverManager<ScalarType,MV,OP> {
123 
124 private:
127  typedef Teuchos::ScalarTraits<ScalarType> SCT;
128  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
129  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
130 
131 public:
132 
134 
135 
142 
163  BlockGmresSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
164  const Teuchos::RCP<Teuchos::ParameterList> &pl );
165 
167  virtual ~BlockGmresSolMgr() {};
168 
170  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
171  return Teuchos::rcp(new BlockGmresSolMgr<ScalarType,MV,OP>);
172  }
174 
176 
177 
180  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
181  return *problem_;
182  }
183 
186  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
187 
190  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
191 
197  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
198  return Teuchos::tuple(timerSolve_);
199  }
200 
211  MagnitudeType achievedTol() const override {
212  return achievedTol_;
213  }
214 
216  int getNumIters() const override {
217  return numIters_;
218  }
219 
223  bool isLOADetected() const override { return loaDetected_; }
224 
226 
228 
229 
231  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; isSTSet_ = false; }
232 
234  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
235 
237  void setDebugStatusTest( const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &debugStatusTest ) override;
238 
240 
242 
243 
247  void reset( const ResetType type ) override { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
249 
251 
252 
270  ReturnType solve() override;
271 
273 
276 
283  void
284  describe (Teuchos::FancyOStream& out,
285  const Teuchos::EVerbosityLevel verbLevel =
286  Teuchos::Describable::verbLevel_default) const override;
287 
289  std::string description () const override;
290 
292 
293 private:
294 
295  // Method for checking current status test against defined linear problem.
296  bool checkStatusTest();
297 
298  // Linear problem.
299  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
300 
301  // Output manager.
302  Teuchos::RCP<OutputManager<ScalarType> > printer_;
303  Teuchos::RCP<std::ostream> outputStream_;
304 
305  // Status test.
306  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugStatusTest_;
307  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
308  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
309  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
310  Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
311  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
312 
313  // Orthogonalization manager.
314  Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
315 
316  // Current parameter list.
317  Teuchos::RCP<Teuchos::ParameterList> params_;
318 
319  // Default solver values.
320  static constexpr int maxRestarts_default_ = 20;
321  static constexpr int maxIters_default_ = 1000;
322  static constexpr bool adaptiveBlockSize_default_ = true;
323  static constexpr bool showMaxResNormOnly_default_ = false;
324  static constexpr bool flexibleGmres_default_ = false;
325  static constexpr bool expResTest_default_ = false;
326  static constexpr int blockSize_default_ = 1;
327  static constexpr int numBlocks_default_ = 300;
328  static constexpr int verbosity_default_ = Belos::Errors;
329  static constexpr int outputStyle_default_ = Belos::General;
330  static constexpr int outputFreq_default_ = -1;
331  static constexpr const char * impResScale_default_ = "Norm of Preconditioned Initial Residual";
332  static constexpr const char * expResScale_default_ = "Norm of Initial Residual";
333  static constexpr const char * label_default_ = "Belos";
334  static constexpr const char * orthoType_default_ = "ICGS";
335 // https://stackoverflow.com/questions/24398102/constexpr-and-initialization-of-a-static-const-void-pointer-with-reinterpret-cas
336 #if defined(_WIN32) && defined(__clang__)
337  static constexpr std::ostream * outputStream_default_ =
338  __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
339 #else
340  static constexpr std::ostream * outputStream_default_ = &std::cout;
341 #endif
342 
343  // Current solver values.
348  std::string orthoType_;
350 
351  // Timers.
352  std::string label_;
353  Teuchos::RCP<Teuchos::Time> timerSolve_;
354 
355  // Internal state variables.
358 };
359 
360 
361 // Empty Constructor
362 template<class ScalarType, class MV, class OP>
364  outputStream_(Teuchos::rcp(outputStream_default_,false)),
365  convtol_(DefaultSolverParameters::convTol),
366  orthoKappa_(DefaultSolverParameters::orthoKappa),
367  achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
368  maxRestarts_(maxRestarts_default_),
369  maxIters_(maxIters_default_),
370  numIters_(0),
371  blockSize_(blockSize_default_),
372  numBlocks_(numBlocks_default_),
373  verbosity_(verbosity_default_),
374  outputStyle_(outputStyle_default_),
375  outputFreq_(outputFreq_default_),
376  adaptiveBlockSize_(adaptiveBlockSize_default_),
377  showMaxResNormOnly_(showMaxResNormOnly_default_),
378  isFlexible_(flexibleGmres_default_),
379  expResTest_(expResTest_default_),
380  orthoType_(orthoType_default_),
381  impResScale_(impResScale_default_),
382  expResScale_(expResScale_default_),
383  label_(label_default_),
384  isSet_(false),
385  isSTSet_(false),
386  loaDetected_(false)
387 {}
388 
389 
390 // Basic Constructor
391 template<class ScalarType, class MV, class OP>
393 BlockGmresSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
394  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
395  problem_(problem),
396  outputStream_(Teuchos::rcp(outputStream_default_,false)),
397  convtol_(DefaultSolverParameters::convTol),
398  orthoKappa_(DefaultSolverParameters::orthoKappa),
399  achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
400  maxRestarts_(maxRestarts_default_),
401  maxIters_(maxIters_default_),
402  numIters_(0),
403  blockSize_(blockSize_default_),
404  numBlocks_(numBlocks_default_),
405  verbosity_(verbosity_default_),
406  outputStyle_(outputStyle_default_),
407  outputFreq_(outputFreq_default_),
408  adaptiveBlockSize_(adaptiveBlockSize_default_),
409  showMaxResNormOnly_(showMaxResNormOnly_default_),
410  isFlexible_(flexibleGmres_default_),
411  expResTest_(expResTest_default_),
412  orthoType_(orthoType_default_),
413  impResScale_(impResScale_default_),
414  expResScale_(expResScale_default_),
415  label_(label_default_),
416  isSet_(false),
417  isSTSet_(false),
418  loaDetected_(false)
419 {
420 
421  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
422 
423  // If the parameter list pointer is null, then set the current parameters to the default parameter list.
424  if ( !is_null(pl) ) {
425  setParameters( pl );
426  }
427 
428 }
429 
430 
431 template<class ScalarType, class MV, class OP>
432 Teuchos::RCP<const Teuchos::ParameterList>
434 {
435  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
436  if (is_null(validPL)) {
437  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
438 
439  // The static_cast is to resolve an issue with older clang versions which
440  // would cause the constexpr to link fail. With c++17 the problem is resolved.
441  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
442  "The relative residual tolerance that needs to be achieved by the\n"
443  "iterative solver in order for the linear system to be declared converged." );
444  pl->set("Maximum Restarts", static_cast<int>(maxRestarts_default_),
445  "The maximum number of restarts allowed for each\n"
446  "set of RHS solved.");
447  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
448  "The maximum number of block iterations allowed for each\n"
449  "set of RHS solved.");
450  pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
451  "The maximum number of blocks allowed in the Krylov subspace\n"
452  "for each set of RHS solved.");
453  pl->set("Block Size", static_cast<int>(blockSize_default_),
454  "The number of vectors in each block. This number times the\n"
455  "number of blocks is the total Krylov subspace dimension.");
456  pl->set("Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
457  "Whether the solver manager should adapt the block size\n"
458  "based on the number of RHS to solve.");
459  pl->set("Verbosity", static_cast<int>(verbosity_default_),
460  "What type(s) of solver information should be outputted\n"
461  "to the output stream.");
462  pl->set("Output Style", static_cast<int>(outputStyle_default_),
463  "What style is used for the solver information outputted\n"
464  "to the output stream.");
465  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
466  "How often convergence information should be outputted\n"
467  "to the output stream.");
468  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
469  "A reference-counted pointer to the output stream where all\n"
470  "solver output is sent.");
471  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
472  "When convergence information is printed, only show the maximum\n"
473  "relative residual norm when the block size is greater than one.");
474  pl->set("Flexible Gmres", static_cast<bool>(flexibleGmres_default_),
475  "Whether the solver manager should use the flexible variant\n"
476  "of GMRES.");
477  pl->set("Explicit Residual Test", static_cast<bool>(expResTest_default_),
478  "Whether the explicitly computed residual should be used in the convergence test.");
479  pl->set("Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
480  "The type of scaling used in the implicit residual convergence test.");
481  pl->set("Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
482  "The type of scaling used in the explicit residual convergence test.");
483  pl->set("Timer Label", static_cast<const char *>(label_default_),
484  "The string to use as a prefix for the timer labels.");
485  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
486  "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
487  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
488  "The constant used by DGKS orthogonalization to determine\n"
489  "whether another step of classical Gram-Schmidt is necessary.");
490  validPL = pl;
491  }
492  return validPL;
493 }
494 
495 
496 template<class ScalarType, class MV, class OP>
497 void BlockGmresSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
498 {
499 
500  // Create the internal parameter list if ones doesn't already exist.
501  if (params_ == Teuchos::null) {
502  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
503  }
504  else {
505  params->validateParameters(*getValidParameters());
506  }
507 
508  // Check for maximum number of restarts
509  if (params->isParameter("Maximum Restarts")) {
510  maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
511 
512  // Update parameter in our list.
513  params_->set("Maximum Restarts", maxRestarts_);
514  }
515 
516  // Check for maximum number of iterations
517  if (params->isParameter("Maximum Iterations")) {
518  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
519 
520  // Update parameter in our list and in status test.
521  params_->set("Maximum Iterations", maxIters_);
522  if (maxIterTest_!=Teuchos::null)
523  maxIterTest_->setMaxIters( maxIters_ );
524  }
525 
526  // Check for blocksize
527  if (params->isParameter("Block Size")) {
528  blockSize_ = params->get("Block Size",blockSize_default_);
529  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
530  "Belos::BlockGmresSolMgr: \"Block Size\" must be strictly positive.");
531 
532  // Update parameter in our list.
533  params_->set("Block Size", blockSize_);
534  }
535 
536  // Check if the blocksize should be adaptive
537  if (params->isParameter("Adaptive Block Size")) {
538  adaptiveBlockSize_ = params->get("Adaptive Block Size",adaptiveBlockSize_default_);
539 
540  // Update parameter in our list.
541  params_->set("Adaptive Block Size", adaptiveBlockSize_);
542  }
543 
544  // Check for the maximum number of blocks.
545  if (params->isParameter("Num Blocks")) {
546  numBlocks_ = params->get("Num Blocks",numBlocks_default_);
547  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
548  "Belos::BlockGmresSolMgr: \"Num Blocks\" must be strictly positive.");
549 
550  // Update parameter in our list.
551  params_->set("Num Blocks", numBlocks_);
552  }
553 
554  // Check to see if the timer label changed.
555  if (params->isParameter("Timer Label")) {
556  std::string tempLabel = params->get("Timer Label", label_default_);
557 
558  // Update parameter in our list, solver timer, and orthogonalization label
559  if (tempLabel != label_) {
560  label_ = tempLabel;
561  params_->set("Timer Label", label_);
562  std::string solveLabel = label_ + ": BlockGmresSolMgr total solve time";
563 #ifdef BELOS_TEUCHOS_TIME_MONITOR
564  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
565 #endif
566  if (ortho_ != Teuchos::null) {
567  ortho_->setLabel( label_ );
568  }
569  }
570  }
571 
572  // Determine whether this solver should be "flexible".
573  if (params->isParameter("Flexible Gmres")) {
574  isFlexible_ = Teuchos::getParameter<bool>(*params,"Flexible Gmres");
575  params_->set("Flexible Gmres", isFlexible_);
576  if (isFlexible_ && expResTest_) {
577  // Use an implicit convergence test if the Gmres solver is flexible
578  isSTSet_ = false;
579  }
580  }
581 
582  // Check for a change in verbosity level
583  if (params->isParameter("Verbosity")) {
584  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
585  verbosity_ = params->get("Verbosity", verbosity_default_);
586  } else {
587  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
588  }
589 
590  // Update parameter in our list.
591  params_->set("Verbosity", verbosity_);
592  if (printer_ != Teuchos::null)
593  printer_->setVerbosity(verbosity_);
594  }
595 
596  // Check for a change in output style
597  if (params->isParameter("Output Style")) {
598  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
599  outputStyle_ = params->get("Output Style", outputStyle_default_);
600  } else {
601  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
602  }
603 
604  // Reconstruct the convergence test if the explicit residual test is not being used.
605  params_->set("Output Style", outputStyle_);
606  if (outputTest_ != Teuchos::null) {
607  isSTSet_ = false;
608  }
609  }
610 
611  // output stream
612  if (params->isParameter("Output Stream")) {
613  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
614 
615  // Update parameter in our list.
616  params_->set("Output Stream", outputStream_);
617  if (printer_ != Teuchos::null)
618  printer_->setOStream( outputStream_ );
619  }
620 
621  // frequency level
622  if (verbosity_ & Belos::StatusTestDetails) {
623  if (params->isParameter("Output Frequency")) {
624  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
625  }
626 
627  // Update parameter in out list and output status test.
628  params_->set("Output Frequency", outputFreq_);
629  if (outputTest_ != Teuchos::null)
630  outputTest_->setOutputFrequency( outputFreq_ );
631  }
632 
633  // Create output manager if we need to.
634  if (printer_ == Teuchos::null) {
635  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
636  }
637 
638  // Check if the orthogonalization changed.
639  bool changedOrthoType = false;
640  if (params->isParameter("Orthogonalization")) {
641  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
642  if (tempOrthoType != orthoType_) {
643  orthoType_ = tempOrthoType;
644  changedOrthoType = true;
645  }
646  }
647  params_->set("Orthogonalization", orthoType_);
648 
649  // Check which orthogonalization constant to use.
650  if (params->isParameter("Orthogonalization Constant")) {
651  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
652  orthoKappa_ = params->get ("Orthogonalization Constant",
653  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
654  }
655  else {
656  orthoKappa_ = params->get ("Orthogonalization Constant",
658  }
659 
660  // Update parameter in our list.
661  params_->set("Orthogonalization Constant",orthoKappa_);
662  if (orthoType_=="DGKS") {
663  if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
664  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
665  }
666  }
667  }
668 
669  // Create orthogonalization manager if we need to.
670  if (ortho_ == Teuchos::null || changedOrthoType) {
672  Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
673  if (orthoType_=="DGKS" && orthoKappa_ > 0) {
674  paramsOrtho->set ("depTol", orthoKappa_ );
675  }
676 
677  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
678  }
679 
680  // Check for convergence tolerance
681  if (params->isParameter("Convergence Tolerance")) {
682  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
683  convtol_ = params->get ("Convergence Tolerance",
684  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
685  }
686  else {
687  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
688  }
689 
690  // Update parameter in our list and residual tests.
691  params_->set("Convergence Tolerance", convtol_);
692  if (impConvTest_ != Teuchos::null)
693  impConvTest_->setTolerance( convtol_ );
694  if (expConvTest_ != Teuchos::null)
695  expConvTest_->setTolerance( convtol_ );
696  }
697 
698  // Check for a change in scaling, if so we need to build new residual tests.
699  if (params->isParameter("Implicit Residual Scaling")) {
700  std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
701 
702  // Only update the scaling if it's different.
703  if (impResScale_ != tempImpResScale) {
704  Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
705  impResScale_ = tempImpResScale;
706 
707  // Update parameter in our list and residual tests
708  params_->set("Implicit Residual Scaling", impResScale_);
709  if (impConvTest_ != Teuchos::null) {
710  try {
711  impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
712  }
713  catch (std::exception& e) {
714  // Make sure the convergence test gets constructed again.
715  isSTSet_ = false;
716  }
717  }
718  }
719  }
720 
721  if (params->isParameter("Explicit Residual Scaling")) {
722  std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
723 
724  // Only update the scaling if it's different.
725  if (expResScale_ != tempExpResScale) {
726  Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
727  expResScale_ = tempExpResScale;
728 
729  // Update parameter in our list and residual tests
730  params_->set("Explicit Residual Scaling", expResScale_);
731  if (expConvTest_ != Teuchos::null) {
732  try {
733  expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
734  }
735  catch (std::exception& e) {
736  // Make sure the convergence test gets constructed again.
737  isSTSet_ = false;
738  }
739  }
740  }
741  }
742 
743  if (params->isParameter("Explicit Residual Test")) {
744  expResTest_ = Teuchos::getParameter<bool>( *params,"Explicit Residual Test" );
745 
746  // Reconstruct the convergence test if the explicit residual test is not being used.
747  params_->set("Explicit Residual Test", expResTest_);
748  if (expConvTest_ == Teuchos::null) {
749  isSTSet_ = false;
750  }
751  }
752 
753  if (params->isParameter("Show Maximum Residual Norm Only")) {
754  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
755 
756  // Update parameter in our list and residual tests
757  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
758  if (impConvTest_ != Teuchos::null)
759  impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
760  if (expConvTest_ != Teuchos::null)
761  expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
762  }
763 
764 
765  // Create the timer if we need to.
766  if (timerSolve_ == Teuchos::null) {
767  std::string solveLabel = label_ + ": BlockGmresSolMgr total solve time";
768 #ifdef BELOS_TEUCHOS_TIME_MONITOR
769  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
770 #endif
771  }
772 
773  // Inform the solver manager that the current parameters were set.
774  isSet_ = true;
775 }
776 
777 // Check the status test versus the defined linear problem
778 template<class ScalarType, class MV, class OP>
780 
781  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
782  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestGenResNorm_t;
783  typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP> StatusTestImpResNorm_t;
784 
785  // Basic test checks maximum iterations and native residual.
786  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
787 
788  // Perform sanity checking for flexible Gmres here.
789  // NOTE: If the user requests that the solver manager use flexible GMRES, but there is no right preconditioner, don't use flexible GMRES.
790  // Throw an error is the user provided a left preconditioner, as that is inconsistent with flexible GMRES.
791  if (isFlexible_ && Teuchos::is_null(problem_->getRightPrec())) {
792  isFlexible_ = false;
793  params_->set("Flexible Gmres", isFlexible_);
794 
795  // If the user specified the preconditioner as a left preconditioner, throw an error.
796  TEUCHOS_TEST_FOR_EXCEPTION( !Teuchos::is_null(problem_->getLeftPrec()),BlockGmresSolMgrLinearProblemFailure,
797  "Belos::BlockGmresSolMgr::solve(): Linear problem has a left preconditioner, not a right preconditioner, which is required for flexible GMRES.");
798  }
799 
800  // If there is a left preconditioner, we create a combined status test that checks the implicit
801  // and then explicit residual norm to see if we have convergence.
802  if (!Teuchos::is_null(problem_->getLeftPrec()) && !isFlexible_) {
803  expResTest_ = true;
804  }
805 
806  if (expResTest_) {
807 
808  // Implicit residual test, using the native residual to determine if convergence was achieved.
809  Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
810  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
811  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
812  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
813  impConvTest_ = tmpImpConvTest;
814 
815  // Explicit residual test once the native residual is below the tolerance
816  Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
817  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
818  tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
819  tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
820  tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
821  expConvTest_ = tmpExpConvTest;
822 
823  // The convergence test is a combination of the "cheap" implicit test and explicit test.
824  convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
825  }
826  else {
827 
828  if (isFlexible_) {
829  // Implicit residual test, using the native residual to determine if convergence was achieved.
830  Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
831  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
832  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
833  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
834  impConvTest_ = tmpImpConvTest;
835  }
836  else {
837  // Implicit residual test, using the native residual to determine if convergence was achieved.
838  // Use test that checks for loss of accuracy.
839  Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
840  Teuchos::rcp( new StatusTestImpResNorm_t( convtol_ ) );
841  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
842  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
843  impConvTest_ = tmpImpConvTest;
844  }
845 
846  // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
847  expConvTest_ = impConvTest_;
848  convTest_ = impConvTest_;
849  }
850 
851  // Create the status test.
852  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
853 
854  // Add debug status test, if one is provided by the user
855  if (nonnull(debugStatusTest_) ) {
856  // Add debug convergence test
857  Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
858  }
859 
860  // Create the status test output class.
861  // This class manages and formats the output from the status test.
862  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
863  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
864 
865  // Set the solver string for the output test
866  std::string solverDesc = " Block Gmres ";
867  if (isFlexible_)
868  solverDesc = "Flexible" + solverDesc;
869  outputTest_->setSolverDesc( solverDesc );
870 
871  // The status test is now set.
872  isSTSet_ = true;
873 
874  return false;
875 }
876 
877 template<class ScalarType, class MV, class OP>
879  const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &debugStatusTest
880  )
881 {
882  debugStatusTest_ = debugStatusTest;
883 }
884 
885 
886 // solve()
887 template<class ScalarType, class MV, class OP>
889 
890  // Set the current parameters if they were not set before.
891  // NOTE: This may occur if the user generated the solver manager with the default constructor and
892  // then didn't set any parameters using setParameters().
893  if (!isSet_) {
894  setParameters(Teuchos::parameterList(*getValidParameters()));
895  }
896 
897  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,BlockGmresSolMgrLinearProblemFailure,
898  "Belos::BlockGmresSolMgr::solve(): Linear problem is not a valid object.");
899 
900  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),BlockGmresSolMgrLinearProblemFailure,
901  "Belos::BlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
902 
903  if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
904  TEUCHOS_TEST_FOR_EXCEPTION( checkStatusTest(),BlockGmresSolMgrLinearProblemFailure,
905  "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
906  }
907 
908  // Create indices for the linear systems to be solved.
909  int startPtr = 0;
910  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
911  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
912 
913  std::vector<int> currIdx;
914  // If an adaptive block size is allowed then only the linear systems that need to be solved are solved.
915  // Otherwise, the index set is generated that informs the linear problem that some linear systems are augmented.
916  if ( adaptiveBlockSize_ ) {
917  blockSize_ = numCurrRHS;
918  currIdx.resize( numCurrRHS );
919  for (int i=0; i<numCurrRHS; ++i)
920  { currIdx[i] = startPtr+i; }
921 
922  }
923  else {
924  currIdx.resize( blockSize_ );
925  for (int i=0; i<numCurrRHS; ++i)
926  { currIdx[i] = startPtr+i; }
927  for (int i=numCurrRHS; i<blockSize_; ++i)
928  { currIdx[i] = -1; }
929  }
930 
931  // Inform the linear problem of the current linear system to solve.
932  problem_->setLSIndex( currIdx );
933 
935  // Parameter list
936  Teuchos::ParameterList plist;
937  plist.set("Block Size",blockSize_);
938 
939  ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
940  if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
941  int tmpNumBlocks = 0;
942  if (blockSize_ == 1)
943  tmpNumBlocks = dim / blockSize_; // Allow for a good breakdown.
944  else
945  tmpNumBlocks = ( dim - blockSize_) / blockSize_; // Allow for restarting.
946  printer_->stream(Warnings) <<
947  "Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!"
948  << std::endl << " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
949  plist.set("Num Blocks",tmpNumBlocks);
950  }
951  else
952  plist.set("Num Blocks",numBlocks_);
953 
954  // Reset the status test.
955  outputTest_->reset();
956  loaDetected_ = false;
957 
958  // Assume convergence is achieved, then let any failed convergence set this to false.
959  bool isConverged = true;
960 
962  // BlockGmres solver
963 
964  Teuchos::RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter;
965 
966  if (isFlexible_)
967  block_gmres_iter = Teuchos::rcp( new BlockFGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
968  else
969  block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
970 
971  // Enter solve() iterations
972  {
973 #ifdef BELOS_TEUCHOS_TIME_MONITOR
974  Teuchos::TimeMonitor slvtimer(*timerSolve_);
975 #endif
976 
977  while ( numRHS2Solve > 0 ) {
978 
979  // Set the current number of blocks and blocksize with the Gmres iteration.
980  if (blockSize_*numBlocks_ > dim) {
981  int tmpNumBlocks = 0;
982  if (blockSize_ == 1)
983  tmpNumBlocks = dim / blockSize_; // Allow for a good breakdown.
984  else
985  tmpNumBlocks = ( dim - blockSize_) / blockSize_; // Allow for restarting.
986  block_gmres_iter->setSize( blockSize_, tmpNumBlocks );
987  }
988  else
989  block_gmres_iter->setSize( blockSize_, numBlocks_ );
990 
991  // Reset the number of iterations.
992  block_gmres_iter->resetNumIters();
993 
994  // Reset the number of calls that the status test output knows about.
995  outputTest_->resetNumCalls();
996 
997  // Create the first block in the current Krylov basis.
998  Teuchos::RCP<MV> V_0;
999  if (isFlexible_) {
1000  // Load the correct residual if the system is augmented
1001  if (currIdx[blockSize_-1] == -1) {
1002  V_0 = MVT::Clone( *(problem_->getInitResVec()), blockSize_ );
1003  problem_->computeCurrResVec( &*V_0 );
1004  }
1005  else {
1006  V_0 = MVT::CloneCopy( *(problem_->getInitResVec()), currIdx );
1007  }
1008  }
1009  else {
1010  // Load the correct residual if the system is augmented
1011  if (currIdx[blockSize_-1] == -1) {
1012  V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1013  problem_->computeCurrPrecResVec( &*V_0 );
1014  }
1015  else {
1016  V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1017  }
1018  }
1019 
1020  // Get a matrix to hold the orthonormalization coefficients.
1021  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 =
1022  Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
1023 
1024  // Orthonormalize the new V_0
1025  int rank = ortho_->normalize( *V_0, z_0 );
1026  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockGmresSolMgrOrthoFailure,
1027  "Belos::BlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1028 
1029  // Set the new state and initialize the solver.
1031  newstate.V = V_0;
1032  newstate.z = z_0;
1033  newstate.curDim = 0;
1034  block_gmres_iter->initializeGmres(newstate);
1035  int numRestarts = 0;
1036 
1037  while(1) {
1038  // tell block_gmres_iter to iterate
1039  try {
1040  block_gmres_iter->iterate();
1041 
1043  //
1044  // check convergence first
1045  //
1047  if ( convTest_->getStatus() == Passed ) {
1048  if ( expConvTest_->getLOADetected() ) {
1049  // we don't have convergence
1050  loaDetected_ = true;
1051  printer_->stream(Warnings) <<
1052  "Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1053  isConverged = false;
1054  }
1055  break; // break from while(1){block_gmres_iter->iterate()}
1056  }
1058  //
1059  // check for maximum iterations
1060  //
1062  else if ( maxIterTest_->getStatus() == Passed ) {
1063  // we don't have convergence
1064  isConverged = false;
1065  break; // break from while(1){block_gmres_iter->iterate()}
1066  }
1068  //
1069  // check for restarting, i.e. the subspace is full
1070  //
1072  else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1073 
1074  if ( numRestarts >= maxRestarts_ ) {
1075  isConverged = false;
1076  break; // break from while(1){block_gmres_iter->iterate()}
1077  }
1078  numRestarts++;
1079 
1080  printer_->stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
1081 
1082  // Update the linear problem.
1083  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1084  if (isFlexible_) {
1085  // Update the solution manually, since the preconditioning doesn't need to be undone.
1086  Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1087  MVT::MvAddMv( 1.0, *curX, 1.0, *update, *curX );
1088  }
1089  else
1090  problem_->updateSolution( update, true );
1091 
1092  // Get the state.
1093  GmresIterationState<ScalarType,MV> oldState = block_gmres_iter->getState();
1094 
1095  // Compute the restart std::vector.
1096  // Get a view of the current Krylov basis.
1097  V_0 = MVT::Clone( *(oldState.V), blockSize_ );
1098  if (isFlexible_)
1099  problem_->computeCurrResVec( &*V_0 );
1100  else
1101  problem_->computeCurrPrecResVec( &*V_0 );
1102 
1103  // Get a view of the first block of the Krylov basis.
1104  z_0 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
1105 
1106  // Orthonormalize the new V_0
1107  rank = ortho_->normalize( *V_0, z_0 );
1108  TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockGmresSolMgrOrthoFailure,
1109  "Belos::BlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after restart.");
1110 
1111  // Set the new state and initialize the solver.
1112  newstate.V = V_0;
1113  newstate.z = z_0;
1114  newstate.curDim = 0;
1115  block_gmres_iter->initializeGmres(newstate);
1116 
1117  } // end of restarting
1118 
1120  //
1121  // we returned from iterate(), but none of our status tests Passed.
1122  // something is wrong, and it is probably our fault.
1123  //
1125 
1126  else {
1127  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1128  "Belos::BlockGmresSolMgr::solve(): Invalid return from BlockGmresIter::iterate().");
1129  }
1130  }
1131  catch (const GmresIterationOrthoFailure &e) {
1132  // If the block size is not one, it's not considered a lucky breakdown.
1133  if (blockSize_ != 1) {
1134  printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
1135  << block_gmres_iter->getNumIters() << std::endl
1136  << e.what() << std::endl;
1137  if (convTest_->getStatus() != Passed)
1138  isConverged = false;
1139  break;
1140  }
1141  else {
1142  // If the block size is one, try to recover the most recent least-squares solution
1143  block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1144 
1145  // Check to see if the most recent least-squares solution yielded convergence.
1146  sTest_->checkStatus( &*block_gmres_iter );
1147  if (convTest_->getStatus() != Passed)
1148  isConverged = false;
1149  break;
1150  }
1151  }
1152  catch (const std::exception &e) {
1153  printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
1154  << block_gmres_iter->getNumIters() << std::endl
1155  << e.what() << std::endl;
1156  throw;
1157  }
1158  }
1159 
1160  // Compute the current solution.
1161  // Update the linear problem.
1162  if (isFlexible_) {
1163  // Update the solution manually, since the preconditioning doesn't need to be undone.
1164  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1165  Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1166  // Update the solution only if there is a valid update from the iteration
1167  if (update != Teuchos::null)
1168  MVT::MvAddMv( 1.0, *curX, 1.0, *update, *curX );
1169  }
1170  else {
1171  // Attempt to get the current solution from the residual status test, if it has one.
1172  if ( !Teuchos::is_null(expConvTest_->getSolution()) ) {
1173  Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1174  Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1175  MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1176  }
1177  else {
1178  Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1179  problem_->updateSolution( update, true );
1180  }
1181  }
1182 
1183  // Inform the linear problem that we are finished with this block linear system.
1184  problem_->setCurrLS();
1185 
1186  // Update indices for the linear systems to be solved.
1187  startPtr += numCurrRHS;
1188  numRHS2Solve -= numCurrRHS;
1189  if ( numRHS2Solve > 0 ) {
1190  numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1191 
1192  if ( adaptiveBlockSize_ ) {
1193  blockSize_ = numCurrRHS;
1194  currIdx.resize( numCurrRHS );
1195  for (int i=0; i<numCurrRHS; ++i)
1196  { currIdx[i] = startPtr+i; }
1197  }
1198  else {
1199  currIdx.resize( blockSize_ );
1200  for (int i=0; i<numCurrRHS; ++i)
1201  { currIdx[i] = startPtr+i; }
1202  for (int i=numCurrRHS; i<blockSize_; ++i)
1203  { currIdx[i] = -1; }
1204  }
1205  // Set the next indices.
1206  problem_->setLSIndex( currIdx );
1207  }
1208  else {
1209  currIdx.resize( numRHS2Solve );
1210  }
1211 
1212  }// while ( numRHS2Solve > 0 )
1213 
1214  }
1215 
1216  // print final summary
1217  sTest_->print( printer_->stream(FinalSummary) );
1218 
1219  // print timing information
1220 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1221  // Calling summarize() can be expensive, so don't call unless the
1222  // user wants to print out timing details. summarize() will do all
1223  // the work even if it's passed a "black hole" output stream.
1224  if (verbosity_ & TimingDetails)
1225  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1226 #endif
1227 
1228  // get iteration information for this solve
1229  numIters_ = maxIterTest_->getNumIters();
1230 
1231  // Save the convergence test value ("achieved tolerance") for this
1232  // solve. This requires a bit more work than for BlockCGSolMgr,
1233  // since for this solver, convTest_ may either be a single residual
1234  // norm test, or a combination of two residual norm tests. In the
1235  // latter case, the master convergence test convTest_ is a SEQ combo
1236  // of the implicit resp. explicit tests. If the implicit test never
1237  // passes, then the explicit test won't ever be executed. This
1238  // manifests as expConvTest_->getTestValue()->size() < 1. We deal
1239  // with this case by using the values returned by
1240  // impConvTest_->getTestValue().
1241  {
1242  // We'll fetch the vector of residual norms one way or the other.
1243  const std::vector<MagnitudeType>* pTestValues = NULL;
1244  if (expResTest_) {
1245  pTestValues = expConvTest_->getTestValue();
1246  if (pTestValues == NULL || pTestValues->size() < 1) {
1247  pTestValues = impConvTest_->getTestValue();
1248  }
1249  }
1250  else {
1251  // Only the implicit residual norm test is being used.
1252  pTestValues = impConvTest_->getTestValue();
1253  }
1254  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1255  "Belos::BlockGmresSolMgr::solve(): The implicit convergence test's "
1256  "getTestValue() method returned NULL. Please report this bug to the "
1257  "Belos developers.");
1258  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1259  "Belos::BlockGmresSolMgr::solve(): The implicit convergence test's "
1260  "getTestValue() method returned a vector of length zero. Please report "
1261  "this bug to the Belos developers.");
1262 
1263  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1264  // achieved tolerances for all vectors in the current solve(), or
1265  // just for the vectors from the last deflation?
1266  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1267  }
1268 
1269  if (!isConverged || loaDetected_) {
1270  return Unconverged; // return from BlockGmresSolMgr::solve()
1271  }
1272  return Converged; // return from BlockGmresSolMgr::solve()
1273 }
1274 
1275 
1276 template<class ScalarType, class MV, class OP>
1278 {
1279  std::ostringstream out;
1280  out << "\"Belos::BlockGmresSolMgr\": {";
1281  if (this->getObjectLabel () != "") {
1282  out << "Label: " << this->getObjectLabel () << ", ";
1283  }
1284  out << "Flexible: " << (isFlexible_ ? "true" : "false")
1285  << ", Num Blocks: " << numBlocks_
1286  << ", Maximum Iterations: " << maxIters_
1287  << ", Maximum Restarts: " << maxRestarts_
1288  << ", Convergence Tolerance: " << convtol_
1289  << "}";
1290  return out.str ();
1291 }
1292 
1293 
1294 template<class ScalarType, class MV, class OP>
1295 void
1297 describe (Teuchos::FancyOStream &out,
1298  const Teuchos::EVerbosityLevel verbLevel) const
1299 {
1300  using Teuchos::TypeNameTraits;
1301  using Teuchos::VERB_DEFAULT;
1302  using Teuchos::VERB_NONE;
1303  using Teuchos::VERB_LOW;
1304  // using Teuchos::VERB_MEDIUM;
1305  // using Teuchos::VERB_HIGH;
1306  // using Teuchos::VERB_EXTREME;
1307  using std::endl;
1308 
1309  // Set default verbosity if applicable.
1310  const Teuchos::EVerbosityLevel vl =
1311  (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1312 
1313  if (vl != VERB_NONE) {
1314  Teuchos::OSTab tab0 (out);
1315 
1316  out << "\"Belos::BlockGmresSolMgr\":" << endl;
1317  Teuchos::OSTab tab1 (out);
1318  out << "Template parameters:" << endl;
1319  {
1320  Teuchos::OSTab tab2 (out);
1321  out << "ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1322  << "MV: " << TypeNameTraits<MV>::name () << endl
1323  << "OP: " << TypeNameTraits<OP>::name () << endl;
1324  }
1325  if (this->getObjectLabel () != "") {
1326  out << "Label: " << this->getObjectLabel () << endl;
1327  }
1328  out << "Flexible: " << (isFlexible_ ? "true" : "false") << endl
1329  << "Num Blocks: " << numBlocks_ << endl
1330  << "Maximum Iterations: " << maxIters_ << endl
1331  << "Maximum Restarts: " << maxRestarts_ << endl
1332  << "Convergence Tolerance: " << convtol_ << endl;
1333  }
1334 }
1335 
1336 
1337 } // end Belos namespace
1338 
1339 #endif /* BELOS_BLOCK_GMRES_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< OutputManager< ScalarType > > printer_
Collection of types and exceptions used within the Belos solvers.
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
static constexpr int maxIters_default_
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.
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debugStatusTest) override
Set a debug status test that will be checked at the same time as the top-level status test...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
static constexpr int blockSize_default_
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:120
Teuchos::RCP< Teuchos::Time > timerSolve_
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
Teuchos::RCP< const MV > V
The current Krylov basis.
Virtual base class for StatusTest that printing status tests.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Get a parameter list containing the valid parameters for this object.
static constexpr int numBlocks_default_
static constexpr bool expResTest_default_
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
BlockGmresSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
Structure to contain pointers to GmresIteration state variables.
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:293
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::ScalarTraits< ScalarType > SCT
Interface to Block GMRES and Flexible GMRES.
A pure virtual class for defining the status tests for the Belos iterative solvers.
virtual ~BlockGmresSolMgr()
Destructor.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > debugStatusTest_
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
A factory class for generating StatusTestOutput objects.
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver&#39;s iterate() routine unti...
Traits class which defines basic operations on multivectors.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params) override
Set the parameters the solver manager should use to solve the linear problem.
Belos::StatusTest for logically combining several status tests.
A Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
Pure virtual base class which describes the basic interface for a solver manager. ...
int curDim
The current dimension of the reduction.
static constexpr const char * impResScale_default_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
BlockGmresSolMgrLinearProblemFailure(const std::string &what_arg)
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
static constexpr bool flexibleGmres_default_
static constexpr int outputFreq_default_
Teuchos::RCP< std::ostream > outputStream_
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
This class implements the block flexible GMRES iteration, where a block Krylov subspace is constructe...
int getNumIters() const override
Get the iteration count for the most recent call to solve().
static constexpr const char * label_default_
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
static constexpr int maxRestarts_default_
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
static constexpr int verbosity_default_
Belos concrete class for performing the block GMRES iteration.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const override
Print the object with the given verbosity level to a FancyOStream.
static constexpr bool adaptiveBlockSize_default_
static constexpr const char * expResScale_default_
static constexpr const char * orthoType_default_
static constexpr int outputStyle_default_
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.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
static constexpr bool showMaxResNormOnly_default_
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
BlockGmresSolMgr()
Empty constructor for BlockGmresSolMgr. This constructor takes no arguments and sets the default valu...
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
static constexpr std::ostream * outputStream_default_
BlockGmresSolMgrOrthoFailure(const std::string &what_arg)
A class for extending the status testing capabilities of Belos via logical combinations.
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > expConvTest_
std::string description() const override
Return a one-line description of this object.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
BlockGmresSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate ortho...
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:283
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
OperatorTraits< ScalarType, MV, OP > OPT
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Teuchos::RCP< Teuchos::ParameterList > params_
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > impConvTest_
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
Belos concrete class for performing the block, flexible GMRES iteration.
MultiVecTraits< ScalarType, MV > MVT
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.