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

Generated for Belos by doxygen 1.8.14