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:
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.
323  static const MagnitudeType polytol_default_;
324  static const MagnitudeType convtol_default_;
325  static const MagnitudeType orthoKappa_default_;
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.
343  MagnitudeType polytol_, convtol_, orthoKappa_;
344  int maxDegree_, maxRestarts_, maxIters_, numIters_;
345  int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_;
346  bool strictConvTol_, showMaxResNormOnly_;
347  std::string orthoType_;
348  std::string impResScale_, expResScale_;
349 
350  // Polynomial storage
351  int poly_dim_;
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.
361  bool isPolyBuilt_;
362  bool isSet_, isSTSet_, expResTest_;
363  bool loaDetected_;
364 
366  mutable Teuchos::RCP<const Teuchos::ParameterList> validPL_;
367 };
368 
369 
370 // Default solver values.
371 template<class ScalarType, class MV, class OP>
372 const typename GmresPolySolMgr<ScalarType,MV,OP>::MagnitudeType
374 
375 template<class ScalarType, class MV, class OP>
376 const typename GmresPolySolMgr<ScalarType,MV,OP>::MagnitudeType
378 
379 template<class ScalarType, class MV, class OP>
380 const typename GmresPolySolMgr<ScalarType,MV,OP>::MagnitudeType
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
std::string description() const
Method to return description of the hybrid block GMRES solver manager.
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...
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Class which manages the output and verbosity of the Belos solvers.
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
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
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.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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< 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.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:200
GmresPolySolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
Pure virtual base class which describes the basic interface for a solver manager. ...
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
A linear system to solve, and its associated information.
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.
Hybrid block GMRES iterative linear solver.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:149
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
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.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
GmresPolySolMgrPolynomialFailure is thrown when their is a problem generating the GMRES polynomial fo...
GmresPolySolMgrPolynomialFailure(const std::string &what_arg)
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 ...
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.
Class which defines basic traits for the operator type.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
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 ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
int getNumIters() const
Get the iteration count for the most recent call to solve().
GmresPolySolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthon...
GmresPolySolMgrOrthoFailure(const std::string &what_arg)
GmresPolySolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.

Generated on Thu Jul 21 2016 14:43:57 for Belos by doxygen 1.8.11