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