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

Generated for Belos by doxygen 1.8.14