Belos  Version of the Day
BelosGmresPolySolMgr.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_GMRES_POLY_SOLMGR_HPP
43 #define BELOS_GMRES_POLY_SOLMGR_HPP
44 
48 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 
52 #include "BelosLinearProblem.hpp"
53 #include "BelosSolverManager.hpp"
54 
55 #include "BelosGmresPolyOp.hpp"
56 #include "BelosGmresIteration.hpp"
57 #include "BelosBlockGmresIter.hpp"
64 #include "BelosStatusTestCombo.hpp"
66 #include "BelosOutputManager.hpp"
67 #include "Teuchos_BLAS.hpp"
68 #include "Teuchos_as.hpp"
69 #ifdef BELOS_TEUCHOS_TIME_MONITOR
70 #include "Teuchos_TimeMonitor.hpp"
71 #endif
72 
73 
74 namespace Belos {
75 
77 
78 
86  GmresPolySolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
87  {}};
88 
96  GmresPolySolMgrPolynomialFailure(const std::string& what_arg) : BelosError(what_arg)
97  {}};
98 
106  GmresPolySolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
107  {}};
108 
154 template<class ScalarType, class MV, class OP>
155 class GmresPolySolMgr : public SolverManager<ScalarType,MV,OP> {
156 private:
160  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
162 
163 public:
164 
166 
167 
173  GmresPolySolMgr();
174 
190 
192  virtual ~GmresPolySolMgr() {};
193 
197  }
199 
201 
202 
205  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
206  return *problem_;
207  }
208 
212 
216 
223  return Teuchos::tuple(timerSolve_, timerPoly_);
224  }
225 
227  int getNumIters() const override {
228  return numIters_;
229  }
230 
234  bool isLOADetected() const override { return loaDetected_; }
235 
237 
239 
240 
242  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; isSTSet_ = false; }
243 
245  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
246 
248 
250 
259  void reset( const ResetType type ) override {
260  if ((type & Belos::Problem) && ! problem_.is_null ()) {
261  problem_->setProblem ();
262  isPolyBuilt_ = false; // Rebuild the GMRES polynomial
263  }
264  }
265 
267 
269 
287  ReturnType solve() override;
288 
290 
293 
295  std::string description() const override;
296 
298 
299 private:
300 
301  // Method for checking current status test against defined linear problem.
302  bool checkStatusTest();
303 
304  // Method for generating GMRES polynomial.
305  bool generatePoly();
306 
307  // Linear problem.
309 
310  // Output manager.
312  Teuchos::RCP<std::ostream> outputStream_;
313 
314  // Status test.
318  Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
320 
321  // Orthogonalization manager.
323 
324  // Current parameter list.
326 
327  // Default solver values.
328  static constexpr int maxDegree_default_ = 25;
329  static constexpr int maxRestarts_default_ = 20;
330  static constexpr int maxIters_default_ = 1000;
331  static constexpr bool strictConvTol_default_ = false;
332  static constexpr bool showMaxResNormOnly_default_ = false;
333  static constexpr int blockSize_default_ = 1;
334  static constexpr int numBlocks_default_ = 300;
335  static constexpr int verbosity_default_ = Belos::Errors;
336  static constexpr int outputStyle_default_ = Belos::General;
337  static constexpr int outputFreq_default_ = -1;
338  static constexpr const char * impResScale_default_ = "Norm of RHS";
339  static constexpr const char * expResScale_default_ = "Norm of RHS";
340  static constexpr const char * label_default_ = "Belos";
341  static constexpr const char * orthoType_default_ = "DGKS";
342  static constexpr std::ostream * outputStream_default_ = &std::cout;
343 
344  // Current solver values.
345  MagnitudeType polytol_, convtol_, orthoKappa_;
346  int maxDegree_, maxRestarts_, maxIters_, numIters_;
347  int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_;
348  bool strictConvTol_, showMaxResNormOnly_;
349  std::string orthoType_;
350  std::string impResScale_, expResScale_;
351 
352  // Polynomial storage
353  int poly_dim_;
357 
358  // Timers.
359  std::string label_;
360  Teuchos::RCP<Teuchos::Time> timerSolve_, timerPoly_;
361 
362  // Internal state variables.
363  bool isPolyBuilt_;
364  bool isSet_, isSTSet_, expResTest_;
365  bool loaDetected_;
366 
369 };
370 
371 
372 template<class ScalarType, class MV, class OP>
374  outputStream_ (outputStream_default_),
375  polytol_ (DefaultSolverParameters::polyTol),
376  convtol_ (DefaultSolverParameters::convTol),
377  orthoKappa_ (DefaultSolverParameters::orthoKappa),
378  maxDegree_ (maxDegree_default_),
379  maxRestarts_ (maxRestarts_default_),
380  maxIters_ (maxIters_default_),
381  numIters_ (0),
382  blockSize_ (blockSize_default_),
383  numBlocks_ (numBlocks_default_),
384  verbosity_ (verbosity_default_),
385  outputStyle_ (outputStyle_default_),
386  outputFreq_ (outputFreq_default_),
387  strictConvTol_ (strictConvTol_default_),
388  showMaxResNormOnly_ (showMaxResNormOnly_default_),
389  orthoType_ (orthoType_default_),
390  impResScale_ (impResScale_default_),
391  expResScale_ (expResScale_default_),
392  poly_dim_ (0),
393  label_ (label_default_),
394  isPolyBuilt_ (false),
395  isSet_ (false),
396  isSTSet_ (false),
397  expResTest_ (false),
398  loaDetected_ (false)
399 {}
400 
401 
402 template<class ScalarType, class MV, class OP>
406  problem_ (problem),
407  outputStream_ (outputStream_default_),
408  polytol_ (DefaultSolverParameters::polyTol),
409  convtol_ (DefaultSolverParameters::convTol),
410  orthoKappa_ (DefaultSolverParameters::orthoKappa),
411  maxDegree_ (maxDegree_default_),
412  maxRestarts_ (maxRestarts_default_),
413  maxIters_ (maxIters_default_),
414  numIters_ (0),
415  blockSize_ (blockSize_default_),
416  numBlocks_ (numBlocks_default_),
417  verbosity_ (verbosity_default_),
418  outputStyle_ (outputStyle_default_),
419  outputFreq_ (outputFreq_default_),
420  strictConvTol_ (strictConvTol_default_),
421  showMaxResNormOnly_ (showMaxResNormOnly_default_),
422  orthoType_ (orthoType_default_),
423  impResScale_ (impResScale_default_),
424  expResScale_ (expResScale_default_),
425  poly_dim_ (0),
426  label_ (label_default_),
427  isPolyBuilt_ (false),
428  isSet_ (false),
429  isSTSet_ (false),
430  expResTest_ (false),
431  loaDetected_ (false)
432 {
434  problem_.is_null (), std::invalid_argument,
435  "Belos::GmresPolySolMgr: The given linear problem is null. "
436  "Please call this constructor with a nonnull LinearProblem argument, "
437  "or call the constructor that does not take a LinearProblem.");
438 
439  // If the input parameter list is null, then the parameters take
440  // default values.
441  if (! pl.is_null ()) {
442  setParameters (pl);
443  }
444 }
445 
446 
447 template<class ScalarType, class MV, class OP>
450 {
451  if (validPL_.is_null ()) {
452  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList ();
453 
454  // The static_cast is to resolve an issue with older clang versions which
455  // would cause the constexpr to link fail. With c++17 the problem is resolved.
456  pl->set("Polynomial Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::polyTol),
457  "The relative residual tolerance that used to construct the GMRES polynomial.");
458  pl->set("Maximum Degree", static_cast<int>(maxDegree_default_),
459  "The maximum degree allowed for any GMRES polynomial.");
460  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
461  "The relative residual tolerance that needs to be achieved by the\n"
462  "iterative solver in order for the linear system to be declared converged." );
463  pl->set("Maximum Restarts", static_cast<int>(maxRestarts_default_),
464  "The maximum number of restarts allowed for each\n"
465  "set of RHS solved.");
466  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
467  "The maximum number of block iterations allowed for each\n"
468  "set of RHS solved.");
469  pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
470  "The maximum number of blocks allowed in the Krylov subspace\n"
471  "for each set of RHS solved.");
472  pl->set("Block Size", static_cast<int>(blockSize_default_),
473  "The number of vectors in each block. This number times the\n"
474  "number of blocks is the total Krylov subspace dimension.");
475  pl->set("Verbosity", static_cast<int>(verbosity_default_),
476  "What type(s) of solver information should be outputted\n"
477  "to the output stream.");
478  pl->set("Output Style", static_cast<int>(outputStyle_default_),
479  "What style is used for the solver information outputted\n"
480  "to the output stream.");
481  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
482  "How often convergence information should be outputted\n"
483  "to the output stream.");
484  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
485  "A reference-counted pointer to the output stream where all\n"
486  "solver output is sent.");
487  pl->set("Strict Convergence", static_cast<bool>(strictConvTol_default_),
488  "After polynomial is applied, whether solver should try to achieve\n"
489  "the relative residual tolerance.");
490  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
491  "When convergence information is printed, only show the maximum\n"
492  "relative residual norm when the block size is greater than one.");
493  pl->set("Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
494  "The type of scaling used in the implicit residual convergence test.");
495  pl->set("Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
496  "The type of scaling used in the explicit residual convergence test.");
497  pl->set("Timer Label", static_cast<const char *>(label_default_),
498  "The string to use as a prefix for the timer labels.");
499  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
500  "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
501  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
502  "The constant used by DGKS orthogonalization to determine\n"
503  "whether another step of classical Gram-Schmidt is necessary.");
504  validPL_ = pl;
505  }
506  return validPL_;
507 }
508 
509 
510 template<class ScalarType, class MV, class OP>
513 {
514  // Create the internal parameter list if ones doesn't already exist.
515  if (params_.is_null ()) {
516  params_ = Teuchos::parameterList (*getValidParameters ());
517  }
518  else {
519  params->validateParameters (*getValidParameters ());
520  }
521 
522  // Check for maximum polynomial degree
523  if (params->isParameter("Maximum Degree")) {
524  maxDegree_ = params->get("Maximum Degree",maxDegree_default_);
525 
526  // Update parameter in our list.
527  params_->set("Maximum Degree", maxDegree_);
528  }
529 
530  // Check for maximum number of restarts
531  if (params->isParameter("Maximum Restarts")) {
532  maxRestarts_ = params->get("Maximum Restarts",maxRestarts_default_);
533 
534  // Update parameter in our list.
535  params_->set("Maximum Restarts", maxRestarts_);
536  }
537 
538  // Check for maximum number of iterations
539  if (params->isParameter("Maximum Iterations")) {
540  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
541 
542  // Update parameter in our list and in status test.
543  params_->set("Maximum Iterations", maxIters_);
544  if (maxIterTest_!=Teuchos::null)
545  maxIterTest_->setMaxIters( maxIters_ );
546  }
547 
548  // Check for blocksize
549  if (params->isParameter("Block Size")) {
550  blockSize_ = params->get("Block Size",blockSize_default_);
551  TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
552  "Belos::GmresPolySolMgr: \"Block Size\" must be strictly positive.");
553 
554  // Update parameter in our list.
555  params_->set("Block Size", blockSize_);
556  }
557 
558  // Check for the maximum number of blocks.
559  if (params->isParameter("Num Blocks")) {
560  numBlocks_ = params->get("Num Blocks",numBlocks_default_);
561  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
562  "Belos::GmresPolySolMgr: \"Num Blocks\" must be strictly positive.");
563 
564  // Update parameter in our list.
565  params_->set("Num Blocks", numBlocks_);
566  }
567 
568  // Check to see if the timer label changed.
569  if (params->isParameter("Timer Label")) {
570  std::string tempLabel = params->get("Timer Label", label_default_);
571 
572  // Update parameter in our list and solver timer
573  if (tempLabel != label_) {
574  label_ = tempLabel;
575  params_->set("Timer Label", label_);
576  std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
577 #ifdef BELOS_TEUCHOS_TIME_MONITOR
578  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
579 #endif
580  std::string polyLabel = label_ + ": GmresPolySolMgr polynomial creation time";
581 #ifdef BELOS_TEUCHOS_TIME_MONITOR
582  timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
583 #endif
584  if (ortho_ != Teuchos::null) {
585  ortho_->setLabel( label_ );
586  }
587  }
588  }
589 
590  // Check if the orthogonalization changed.
591  if (params->isParameter("Orthogonalization")) {
592  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
593  TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS",
594  std::invalid_argument,
595  "Belos::GmresPolySolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
596  if (tempOrthoType != orthoType_) {
597  params_->set("Orthogonalization", orthoType_);
598  orthoType_ = tempOrthoType;
599  // Create orthogonalization manager
600  if (orthoType_=="DGKS") {
601  if (orthoKappa_ <= 0) {
602  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
603  }
604  else {
605  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
606  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
607  }
608  }
609  else if (orthoType_=="ICGS") {
610  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
611  }
612  else if (orthoType_=="IMGS") {
613  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
614  }
615  }
616  }
617 
618  // Check which orthogonalization constant to use.
619  if (params->isParameter("Orthogonalization Constant")) {
620  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
621  orthoKappa_ = params->get ("Orthogonalization Constant",
622  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
623  }
624  else {
625  orthoKappa_ = params->get ("Orthogonalization Constant",
627  }
628 
629  // Update parameter in our list.
630  params_->set("Orthogonalization Constant",orthoKappa_);
631  if (orthoType_=="DGKS") {
632  if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
633  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
634  }
635  }
636  }
637 
638  // Check for a change in verbosity level
639  if (params->isParameter("Verbosity")) {
640  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
641  verbosity_ = params->get("Verbosity", verbosity_default_);
642  } else {
643  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
644  }
645 
646  // Update parameter in our list.
647  params_->set("Verbosity", verbosity_);
648  if (printer_ != Teuchos::null)
649  printer_->setVerbosity(verbosity_);
650  }
651 
652  // Check for a change in output style
653  if (params->isParameter("Output Style")) {
654  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
655  outputStyle_ = params->get("Output Style", outputStyle_default_);
656  } else {
657  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
658  }
659 
660  // Reconstruct the convergence test if the explicit residual test is not being used.
661  params_->set("Output Style", outputStyle_);
662  if (outputTest_ != Teuchos::null) {
663  isSTSet_ = false;
664  }
665  }
666 
667  // output stream
668  if (params->isParameter("Output Stream")) {
669  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
670 
671  // Update parameter in our list.
672  params_->set("Output Stream", outputStream_);
673  if (printer_ != Teuchos::null)
674  printer_->setOStream( outputStream_ );
675  }
676 
677  // frequency level
678  if (verbosity_ & Belos::StatusTestDetails) {
679  if (params->isParameter("Output Frequency")) {
680  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
681  }
682 
683  // Update parameter in out list and output status test.
684  params_->set("Output Frequency", outputFreq_);
685  if (outputTest_ != Teuchos::null)
686  outputTest_->setOutputFrequency( outputFreq_ );
687  }
688 
689  // Create output manager if we need to.
690  if (printer_ == Teuchos::null) {
691  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
692  }
693 
694  // Convergence
695  //typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t; // unused
696  //typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t; // unused
697 
698  // Check for polynomial convergence tolerance
699  if (params->isParameter("Polynomial Tolerance")) {
700  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
701  polytol_ = params->get ("Polynomial Tolerance",
702  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
703  }
704  else {
705  polytol_ = params->get ("Polynomial Tolerance", DefaultSolverParameters::convTol);
706  }
707 
708  // Update parameter in our list and residual tests.
709  params_->set("Polynomial Tolerance", polytol_);
710  }
711 
712  // Check for convergence tolerance
713  if (params->isParameter("Convergence Tolerance")) {
714  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
715  convtol_ = params->get ("Convergence Tolerance",
716  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
717  }
718  else {
719  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
720  }
721 
722  // Update parameter in our list and residual tests.
723  params_->set("Convergence Tolerance", convtol_);
724  if (impConvTest_ != Teuchos::null)
725  impConvTest_->setTolerance( convtol_ );
726  if (expConvTest_ != Teuchos::null)
727  expConvTest_->setTolerance( convtol_ );
728  }
729 
730  // Check if user requires solver to reach convergence tolerance
731  if (params->isParameter("Strict Convergence")) {
732  strictConvTol_ = params->get("Strict Convergence",strictConvTol_default_);
733 
734  // Update parameter in our list and residual tests
735  params_->set("Strict Convergence", strictConvTol_);
736  }
737 
738  // Check for a change in scaling, if so we need to build new residual tests.
739  if (params->isParameter("Implicit Residual Scaling")) {
740  std::string tempImpResScale = Teuchos::getParameter<std::string>( *params, "Implicit Residual Scaling" );
741 
742  // Only update the scaling if it's different.
743  if (impResScale_ != tempImpResScale) {
744  Belos::ScaleType impResScaleType = convertStringToScaleType( tempImpResScale );
745  impResScale_ = tempImpResScale;
746 
747  // Update parameter in our list and residual tests
748  params_->set("Implicit Residual Scaling", impResScale_);
749  if (impConvTest_ != Teuchos::null) {
750  try {
751  impConvTest_->defineScaleForm( impResScaleType, Belos::TwoNorm );
752  }
753  catch (std::exception& e) {
754  // Make sure the convergence test gets constructed again.
755  isSTSet_ = false;
756  }
757  }
758  }
759  }
760 
761  if (params->isParameter("Explicit Residual Scaling")) {
762  std::string tempExpResScale = Teuchos::getParameter<std::string>( *params, "Explicit Residual Scaling" );
763 
764  // Only update the scaling if it's different.
765  if (expResScale_ != tempExpResScale) {
766  Belos::ScaleType expResScaleType = convertStringToScaleType( tempExpResScale );
767  expResScale_ = tempExpResScale;
768 
769  // Update parameter in our list and residual tests
770  params_->set("Explicit Residual Scaling", expResScale_);
771  if (expConvTest_ != Teuchos::null) {
772  try {
773  expConvTest_->defineScaleForm( expResScaleType, Belos::TwoNorm );
774  }
775  catch (std::exception& e) {
776  // Make sure the convergence test gets constructed again.
777  isSTSet_ = false;
778  }
779  }
780  }
781  }
782 
783 
784  if (params->isParameter("Show Maximum Residual Norm Only")) {
785  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
786 
787  // Update parameter in our list and residual tests
788  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
789  if (impConvTest_ != Teuchos::null)
790  impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
791  if (expConvTest_ != Teuchos::null)
792  expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
793  }
794 
795  // Create orthogonalization manager if we need to.
796  if (ortho_ == Teuchos::null) {
797  params_->set("Orthogonalization", orthoType_);
798  if (orthoType_=="DGKS") {
799  if (orthoKappa_ <= 0) {
800  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
801  }
802  else {
803  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
804  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
805  }
806  }
807  else if (orthoType_=="ICGS") {
808  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
809  }
810  else if (orthoType_=="IMGS") {
811  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
812  }
813  else {
814  TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
815  "Belos::GmresPolySolMgr(): Invalid orthogonalization type.");
816  }
817  }
818 
819  // Create the timers if we need to.
820  if (timerSolve_ == Teuchos::null) {
821  std::string solveLabel = label_ + ": GmresPolySolMgr total solve time";
822 #ifdef BELOS_TEUCHOS_TIME_MONITOR
823  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
824 #endif
825  }
826 
827  if (timerPoly_ == Teuchos::null) {
828  std::string polyLabel = label_ + ": GmresPolySolMgr polynomial creation time";
829 #ifdef BELOS_TEUCHOS_TIME_MONITOR
830  timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
831 #endif
832  }
833 
834  // Inform the solver manager that the current parameters were set.
835  isSet_ = true;
836 }
837 
838 // Check the status test versus the defined linear problem
839 template<class ScalarType, class MV, class OP>
841 
842  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
843  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestGenResNorm_t;
844  typedef Belos::StatusTestImpResNorm<ScalarType,MV,OP> StatusTestImpResNorm_t;
845 
846  // Basic test checks maximum iterations and native residual.
847  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
848 
849  // If there is a left preconditioner, we create a combined status test that checks the implicit
850  // and then explicit residual norm to see if we have convergence.
851  if (!Teuchos::is_null(problem_->getLeftPrec())) {
852  expResTest_ = true;
853  }
854 
855  if (expResTest_) {
856 
857  // Implicit residual test, using the native residual to determine if convergence was achieved.
858  Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
859  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
860  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
861  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
862  impConvTest_ = tmpImpConvTest;
863 
864  // Explicit residual test once the native residual is below the tolerance
865  Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
866  Teuchos::rcp( new StatusTestGenResNorm_t( convtol_ ) );
867  tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit, Belos::TwoNorm );
868  tmpExpConvTest->defineScaleForm( convertStringToScaleType(expResScale_), Belos::TwoNorm );
869  tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
870  expConvTest_ = tmpExpConvTest;
871 
872  // The convergence test is a combination of the "cheap" implicit test and explicit test.
873  convTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
874  }
875  else {
876 
877  // Implicit residual test, using the native residual to determine if convergence was achieved.
878  // Use test that checks for loss of accuracy.
879  Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
880  Teuchos::rcp( new StatusTestImpResNorm_t( convtol_ ) );
881  tmpImpConvTest->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
882  tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
883  impConvTest_ = tmpImpConvTest;
884 
885  // Set the explicit and total convergence test to this implicit test that checks for accuracy loss.
886  expConvTest_ = impConvTest_;
887  convTest_ = impConvTest_;
888  }
889 
890  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
891 
892  // Create the status test output class.
893  // This class manages and formats the output from the status test.
894  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
895  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
896 
897  // Set the solver string for the output test
898  std::string solverDesc = " Gmres Polynomial ";
899  outputTest_->setSolverDesc( solverDesc );
900 
901 
902  // The status test is now set.
903  isSTSet_ = true;
904 
905  return false;
906 }
907 
908 template<class ScalarType, class MV, class OP>
909 bool GmresPolySolMgr<ScalarType,MV,OP>::generatePoly()
910 {
911  // Create a copy of the linear problem that has a zero initial guess and random RHS.
912  Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
913  Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
914  MVT::MvInit( *newX, STS::zero() );
915  MVT::MvRandom( *newB );
917  Teuchos::rcp( new LinearProblem<ScalarType,MV,OP>( problem_->getOperator(), newX, newB ) );
918  newProblem->setLeftPrec( problem_->getLeftPrec() );
919  newProblem->setRightPrec( problem_->getRightPrec() );
920  newProblem->setLabel("Belos GMRES Poly Generation");
921  newProblem->setProblem();
922  std::vector<int> idx(1,0); // Must set the index to be the first vector (0)!
923  newProblem->setLSIndex( idx );
924 
925  // Create a parameter list for the GMRES iteration.
926  Teuchos::ParameterList polyList;
927 
928  // Tell the block solver that the block size is one.
929  polyList.set("Num Blocks",maxDegree_);
930  polyList.set("Block Size",1);
931  polyList.set("Keep Hessenberg", true);
932 
933  // Create a simple status test that either reaches the relative residual tolerance or maximum polynomial size.
935  Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxDegree_ ) );
936 
937  // Implicit residual test, using the native residual to determine if convergence was achieved.
939  Teuchos::rcp( new StatusTestGenResNorm<ScalarType,MV,OP>( polytol_ ) );
940  convTst->defineScaleForm( convertStringToScaleType(impResScale_), Belos::TwoNorm );
941 
942  // Convergence test that stops the iteration when either are satisfied.
944  Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, maxItrTst, convTst ) );
945 
946  // Create Gmres iteration object to perform one cycle of Gmres.
948  gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(newProblem,printer_,polyTest,ortho_,polyList) );
949 
950  // Create the first block in the current Krylov basis (residual).
951  Teuchos::RCP<MV> V_0 = MVT::Clone( *(newProblem->getRHS()), 1 );
952  newProblem->computeCurrPrecResVec( &*V_0 );
953 
954  // Get a matrix to hold the orthonormalization coefficients.
956 
957  // Orthonormalize the new V_0
958  int rank = ortho_->normalize( *V_0, poly_r0_ );
959  TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GmresPolySolMgrOrthoFailure,
960  "Belos::GmresPolySolMgr::generatePoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
961 
962  // Set the new state and initialize the solver.
963  GmresIterationState<ScalarType,MV> newstate;
964  newstate.V = V_0;
965  newstate.z = poly_r0_;
966  newstate.curDim = 0;
967  gmres_iter->initializeGmres(newstate);
968 
969  // Perform Gmres iteration
970  //bool polyConverged = false; // mfh 30 Sep 2017: unused
971  try {
972  gmres_iter->iterate();
973 
974  // Check convergence first
975  if ( convTst->getStatus() == Passed ) {
976  // we have convergence
977  //polyConverged = true; // mfh 30 Sep 2017: unused
978  }
979  }
980  catch (GmresIterationOrthoFailure e) {
981  // Try to recover the most recent least-squares solution
982  gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
983 
984  // Check to see if the most recent least-squares solution yielded convergence.
985  polyTest->checkStatus( &*gmres_iter );
986  if (convTst->getStatus() == Passed) {
987  //polyConverged = true; // mfh 30 Sep 2017: unused
988  }
989  }
990  catch (std::exception e) {
991  using std::endl;
992  printer_->stream(Errors) << "Error! Caught exception in "
993  "BlockGmresIter::iterate() at iteration " << gmres_iter->getNumIters()
994  << endl << e.what () << endl;
995  throw;
996  }
997 
998  // FIXME (mfh 27 Apr 2013) Why aren't we using polyConverged after
999  // setting it? I'm tired of the compile warning so I'll silence it
1000  // here, but I'm curious why we aren't using the variable.
1001  //(void) polyConverged; // mfh 30 Sep 2017: unused
1002 
1003  // Get the solution for this polynomial, use in comparison below
1004  Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
1005 
1006  // Record polynomial info, get current GMRES state
1007  GmresIterationState<ScalarType,MV> gmresState = gmres_iter->getState();
1008 
1009  // If the polynomial has no dimension, the tolerance is too low, return false
1010  poly_dim_ = gmresState.curDim;
1011  if (poly_dim_ == 0) {
1012  return false;
1013  }
1014  //
1015  // Make a view and then copy the RHS of the least squares problem.
1016  //
1017  poly_y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *gmresState.z, poly_dim_, 1 ) );
1018  poly_H_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( *gmresState.H ) );
1019  //
1020  // Solve the least squares problem.
1021  //
1022  const ScalarType one = STS::one ();
1025  Teuchos::NON_UNIT_DIAG, poly_dim_, 1, one,
1026  gmresState.R->values(), gmresState.R->stride(),
1027  poly_y_->values(), poly_y_->stride() );
1028  //
1029  // Generate the polynomial operator
1030  //
1031  poly_Op_ = Teuchos::rcp(
1032  new Belos::GmresPolyOp<ScalarType,MV,OP>( problem_, poly_H_, poly_y_, poly_r0_ ) );
1033 
1034  return true;
1035 }
1036 
1037 
1038 template<class ScalarType, class MV, class OP>
1040 {
1041  using Teuchos::RCP;
1042  using Teuchos::rcp;
1044 
1045  // Set the current parameters if they were not set before. NOTE:
1046  // This may occur if the user generated the solver manager with the
1047  // default constructor and then didn't set any parameters using
1048  // setParameters().
1049  if (! isSet_) {
1050  setParameters (Teuchos::parameterList (*getValidParameters ()));
1051  }
1052 
1054  problem_.is_null (), GmresPolySolMgrLinearProblemFailure,
1055  "Belos::GmresPolySolMgr::solve: The linear problem has not been set yet, "
1056  "or was set to null. Please call setProblem() with a nonnull input before "
1057  "calling solve().");
1058 
1060  ! problem_->isProblemSet (), GmresPolySolMgrLinearProblemFailure,
1061  "Belos::GmresPolySolMgr::solve: The linear problem is not ready. Please "
1062  "call setProblem() on the LinearProblem object before calling solve().");
1063 
1064  if (! isSTSet_ || (! expResTest_ && ! problem_->getLeftPrec ().is_null ())) {
1065  // checkStatusTest() shouldn't have side effects, but it's still
1066  // not such a good idea to put possibly side-effect-y function
1067  // calls in a macro invocation. It could be expensive if the
1068  // macro evaluates the argument more than once, for example.
1069  const bool check = checkStatusTest ();
1072  "Belos::GmresPolySolMgr::solve: Linear problem and requested status "
1073  "tests are incompatible.");
1074  }
1075 
1076  // If the GMRES polynomial has not been constructed for this
1077  // (nmatrix, preconditioner) pair, generate it.
1078  if (! isPolyBuilt_) {
1079 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1080  Teuchos::TimeMonitor slvtimer(*timerPoly_);
1081 #endif
1082  isPolyBuilt_ = generatePoly();
1083  TEUCHOS_TEST_FOR_EXCEPTION( !isPolyBuilt_ && poly_dim_==0, GmresPolySolMgrPolynomialFailure,
1084  "Belos::GmresPolySolMgr::generatePoly(): Failed to generate a non-trivial polynomial, reduce polynomial tolerance.");
1086  "Belos::GmresPolySolMgr::generatePoly(): Failed to generate polynomial that satisfied requirements.");
1087  }
1088 
1089  // Assume convergence is achieved if user does not require strict convergence.
1090  bool isConverged = true;
1091 
1092  // Solve the linear system using the polynomial
1093  {
1094 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1095  Teuchos::TimeMonitor slvtimer(*timerSolve_);
1096 #endif
1097 
1098  // Apply the polynomial to the current linear system
1099  poly_Op_->Apply( *problem_->getRHS(), *problem_->getLHS() );
1100 
1101  // Reset the problem to acknowledge the updated solution
1102  problem_->setProblem ();
1103 
1104  // If we have to strictly adhere to the requested convergence tolerance, set up a standard GMRES solver.
1105  if (strictConvTol_) {
1106 
1107  // Create indices for the linear systems to be solved.
1108  int startPtr = 0;
1109  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1110  int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1111 
1112 
1113  // If an adaptive block size is allowed then only the linear
1114  // systems that need to be solved are solved. Otherwise, the
1115  // index set is generated that informs the linear problem that
1116  // some linear systems are augmented.
1117 
1118  std::vector<int> currIdx (blockSize_);
1119  for (int i = 0; i < numCurrRHS; ++i) {
1120  currIdx[i] = startPtr+i;
1121  }
1122  for (int i = numCurrRHS; i < blockSize_; ++i) {
1123  currIdx[i] = -1;
1124  }
1125 
1126  // Inform the linear problem of the current linear system to solve.
1127  problem_->setLSIndex (currIdx);
1128 
1130  // Parameter list
1131  Teuchos::ParameterList plist;
1132  plist.set ("Block Size", blockSize_);
1133 
1134  ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1135  if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
1136  ptrdiff_t tmpNumBlocks = 0;
1137  if (blockSize_ == 1) {
1138  tmpNumBlocks = dim / blockSize_; // Allow for a good breakdown.
1139  } else {
1140  tmpNumBlocks = ( dim - blockSize_) / blockSize_; // Allow for restarting.
1141  }
1142  printer_->stream(Warnings)
1143  << "Warning! Requested Krylov subspace dimension is larger than "
1144  << "operator dimension!" << std::endl << "The maximum number of "
1145  << "blocks allowed for the Krylov subspace will be adjusted to "
1146  << tmpNumBlocks << std::endl;
1147  plist.set ("Num Blocks", Teuchos::asSafe<int> (tmpNumBlocks));
1148  }
1149  else {
1150  plist.set ("Num Blocks", numBlocks_);
1151  }
1152 
1153  // Reset the status test.
1154  outputTest_->reset ();
1155  loaDetected_ = false;
1156 
1157  // Assume convergence is achieved, then let any failed
1158  // convergence set this to false.
1159  isConverged = true;
1160 
1161  //
1162  // Solve using BlockGmres
1163  //
1164 
1165  RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter =
1166  rcp (new BlockGmresIter<ScalarType,MV,OP> (problem_, printer_,
1167  outputTest_, ortho_, plist));
1168 
1169  // Enter solve() iterations
1170  while (numRHS2Solve > 0) {
1171  // Set the current number of blocks and block size with the
1172  // Gmres iteration.
1173  if (blockSize_*numBlocks_ > dim) {
1174  int tmpNumBlocks = 0;
1175  if (blockSize_ == 1) {
1176  // Leave space for happy breakdown.
1177  tmpNumBlocks = dim / blockSize_;
1178  } else {
1179  // Leave space for restarting.
1180  tmpNumBlocks = (dim - blockSize_) / blockSize_;
1181  }
1182  block_gmres_iter->setSize (blockSize_, tmpNumBlocks);
1183  }
1184  else {
1185  block_gmres_iter->setSize (blockSize_, numBlocks_);
1186  }
1187 
1188  // Reset the number of iterations.
1189  block_gmres_iter->resetNumIters ();
1190 
1191  // Reset the number of calls that the status test output knows about.
1192  outputTest_->resetNumCalls ();
1193 
1194  // Create the first block in the current Krylov basis.
1195  RCP<MV> V_0 = MVT::CloneCopy (*(problem_->getInitPrecResVec ()), currIdx);
1196 
1197  // Get a matrix to hold the orthonormalization coefficients.
1198  RCP<SDM> z_0 = rcp (new SDM (blockSize_, blockSize_));
1199 
1200  // Orthonormalize the new V_0
1201  int rank = ortho_->normalize (*V_0, z_0);
1203  rank != blockSize_, GmresPolySolMgrOrthoFailure,
1204  "Belos::GmresPolySolMgr::solve: Failed to compute initial block of "
1205  "orthonormal vectors.");
1206 
1207  // Set the new state and initialize the solver.
1209  newstate.V = V_0;
1210  newstate.z = z_0;
1211  newstate.curDim = 0;
1212  block_gmres_iter->initializeGmres(newstate);
1213  int numRestarts = 0;
1214 
1215  while(1) {
1216  try {
1217  block_gmres_iter->iterate();
1218 
1219  // check convergence first
1220  if ( convTest_->getStatus() == Passed ) {
1221  if ( expConvTest_->getLOADetected() ) {
1222  // we don't have convergence
1223  loaDetected_ = true;
1224  isConverged = false;
1225  }
1226  break; // break from while(1){block_gmres_iter->iterate()}
1227  }
1228 
1229  // check for maximum iterations
1230  else if ( maxIterTest_->getStatus() == Passed ) {
1231  // we don't have convergence
1232  isConverged = false;
1233  break; // break from while(1){block_gmres_iter->iterate()}
1234  }
1235  // check for restarting, i.e. the subspace is full
1236  else if (block_gmres_iter->getCurSubspaceDim () ==
1237  block_gmres_iter->getMaxSubspaceDim ()) {
1238  if (numRestarts >= maxRestarts_) {
1239  isConverged = false;
1240  break; // break from while(1){block_gmres_iter->iterate()}
1241  }
1242  numRestarts++;
1243 
1244  printer_->stream(Debug)
1245  << " Performing restart number " << numRestarts << " of "
1246  << maxRestarts_ << std::endl << std::endl;
1247 
1248  // Update the linear problem.
1249  RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1250  problem_->updateSolution (update, true);
1251 
1252  // Get the state.
1253  GmresIterationState<ScalarType,MV> oldState = block_gmres_iter->getState();
1254 
1255  // Compute the restart vector.
1256  // Get a view of the current Krylov basis.
1257  //
1258  // We call this VV_0 to avoid shadowing the previously
1259  // defined V_0 above.
1260  RCP<MV> VV_0 = MVT::Clone (*(oldState.V), blockSize_);
1261  problem_->computeCurrPrecResVec (&*VV_0);
1262 
1263  // Get a view of the first block of the Krylov basis.
1264  //
1265  // We call this zz_0 to avoid shadowing the previously
1266  // defined z_0 above.
1267  RCP<SDM> zz_0 = rcp (new SDM (blockSize_, blockSize_));
1268 
1269  // Orthonormalize the new VV_0
1270  const int theRank = ortho_->normalize( *VV_0, zz_0 );
1272  theRank != blockSize_, GmresPolySolMgrOrthoFailure,
1273  "Belos::GmresPolySolMgr::solve: Failed to compute initial "
1274  "block of orthonormal vectors after restart.");
1275 
1276  // Set the new state and initialize the solver.
1278  theNewState.V = VV_0;
1279  theNewState.z = zz_0;
1280  theNewState.curDim = 0;
1281  block_gmres_iter->initializeGmres (theNewState);
1282  } // end of restarting
1283  //
1284  // We returned from iterate(), but none of our status
1285  // tests Passed. Something is wrong, and it is probably
1286  // our fault.
1287  //
1288  else {
1290  true, std::logic_error,
1291  "Belos::GmresPolySolMgr::solve: Invalid return from "
1292  "BlockGmresIter::iterate(). Please report this bug "
1293  "to the Belos developers.");
1294  }
1295  }
1296  catch (const GmresIterationOrthoFailure& e) {
1297  // If the block size is not one, it's not considered a lucky breakdown.
1298  if (blockSize_ != 1) {
1299  printer_->stream(Errors)
1300  << "Error! Caught std::exception in BlockGmresIter::iterate() "
1301  << "at iteration " << block_gmres_iter->getNumIters()
1302  << std::endl << e.what() << std::endl;
1303  if (convTest_->getStatus() != Passed) {
1304  isConverged = false;
1305  }
1306  break;
1307  }
1308  else {
1309  // If the block size is one, try to recover the most
1310  // recent least-squares solution
1311  block_gmres_iter->updateLSQR (block_gmres_iter->getCurSubspaceDim ());
1312 
1313  // Check to see if the most recent least-squares
1314  // solution yielded convergence.
1315  sTest_->checkStatus (&*block_gmres_iter);
1316  if (convTest_->getStatus() != Passed) {
1317  isConverged = false;
1318  }
1319  break;
1320  }
1321  }
1322  catch (const std::exception &e) {
1323  printer_->stream(Errors)
1324  << "Error! Caught std::exception in BlockGmresIter::iterate() "
1325  << "at iteration " << block_gmres_iter->getNumIters() << std::endl
1326  << e.what() << std::endl;
1327  throw;
1328  }
1329  }
1330 
1331  // Compute the current solution. Update the linear problem.
1332  // Attempt to get the current solution from the residual
1333  // status test, if it has one.
1334  if (! Teuchos::is_null (expConvTest_->getSolution ()) ) {
1335  RCP<MV> newX = expConvTest_->getSolution ();
1336  RCP<MV> curX = problem_->getCurrLHSVec ();
1337  MVT::MvAddMv (STS::zero (), *newX, STS::one (), *newX, *curX);
1338  }
1339  else {
1340  RCP<MV> update = block_gmres_iter->getCurrentUpdate ();
1341  problem_->updateSolution (update, true);
1342  }
1343 
1344  // Inform the linear problem that we are finished with this block linear system.
1345  problem_->setCurrLS ();
1346 
1347  // Update indices for the linear systems to be solved.
1348  startPtr += numCurrRHS;
1349  numRHS2Solve -= numCurrRHS;
1350  if (numRHS2Solve > 0) {
1351  numCurrRHS = (numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1352 
1353  currIdx.resize (blockSize_);
1354  for (int i=0; i<numCurrRHS; ++i) {
1355  currIdx[i] = startPtr+i;
1356  }
1357  for (int i=numCurrRHS; i<blockSize_; ++i) {
1358  currIdx[i] = -1;
1359  }
1360 
1361  // Set the next indices.
1362  problem_->setLSIndex( currIdx );
1363  }
1364  else {
1365  currIdx.resize( numRHS2Solve );
1366  }
1367 
1368  }// while ( numRHS2Solve > 0 )
1369 
1370  // print final summary
1371  sTest_->print( printer_->stream(FinalSummary) );
1372 
1373  } // if (strictConvTol_)
1374  } // timing block
1375 
1376  // print timing information
1377 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1378  // Calling summarize() can be expensive, so don't call unless the
1379  // user wants to print out timing details. summarize() will do all
1380  // the work even if it's passed a "black hole" output stream.
1381  if (verbosity_ & TimingDetails)
1382  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1383 #endif
1384 
1385  if (!isConverged || loaDetected_) {
1386  return Unconverged; // return from GmresPolySolMgr::solve()
1387  }
1388  return Converged; // return from GmresPolySolMgr::solve()
1389 }
1390 
1391 
1392 template<class ScalarType, class MV, class OP>
1394 {
1395  std::ostringstream out;
1396 
1397  out << "\"Belos::GmresPolySolMgr\": {"
1398  << "ScalarType: " << Teuchos::TypeNameTraits<ScalarType>::name ()
1399  << ", Ortho Type: " << orthoType_
1400  << ", Block Size: " << blockSize_
1401  << ", Num Blocks: " << numBlocks_
1402  << ", Max Restarts: " << maxRestarts_;
1403  out << "}";
1404  return out.str ();
1405 }
1406 
1407 } // namespace Belos
1408 
1409 #endif // BELOS_GMRES_POLY_SOLMGR_HPP
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
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.
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver&#39;s iterate() routine unti...
virtual ~GmresPolySolMgr()
Destructor.
Defines the GMRES polynomial operator hybrid-GMRES iterative linear solver.
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:119
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
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.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
Structure to contain pointers to GmresIteration state variables.
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
bool is_null() const
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:205
GmresPolySolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
static constexpr double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:294
int curDim
The current dimension of the reduction.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
A linear system to solve, and its associated information.
static constexpr double polyTol
Relative residual tolerance for matrix polynomial construction.
Definition: BelosTypes.hpp:297
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Class which describes the linear problem to be solved by the iterative solver.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
bool isType(const std::string &name) const
Hybrid block GMRES iterative linear solver.
int getNumIters() const override
Get the iteration count for the most recent call to solve().
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&#39;s class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
void reset(const ResetType type) override
Reset the solver.
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 ...
std::string description() const override
Method to return description of the hybrid block GMRES solver manager.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
GmresPolySolMgrPolynomialFailure is thrown when their is a problem generating the GMRES polynomial fo...
GmresPolySolMgrPolynomialFailure(const std::string &what_arg)
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params) override
Set the parameters the solver manager should use to solve the linear problem.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Get a parameter list containing the valid parameters for this object.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:286
GmresPolySolMgr()
Empty constructor for GmresPolySolMgr. This constructor takes no arguments and sets the default value...
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< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
GmresPolySolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthon...
GmresPolySolMgrOrthoFailure(const std::string &what_arg)
GmresPolySolMgrLinearProblemFailure(const std::string &what_arg)
static std::string name()

Generated for Belos by doxygen 1.8.14