Belos Package Browser (Single Doxygen Collection)  Development
BelosPCPGSolMgr.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_PCPG_SOLMGR_HPP
43 #define BELOS_PCPG_SOLMGR_HPP
44 
48 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 #include "BelosLinearProblem.hpp"
52 #include "BelosSolverManager.hpp"
53 
54 #include "BelosPCPGIter.hpp"
55 
59 #include "BelosStatusTestCombo.hpp"
61 #include "BelosOutputManager.hpp"
62 #include "Teuchos_LAPACK.hpp"
63 #ifdef BELOS_TEUCHOS_TIME_MONITOR
64 # include "Teuchos_TimeMonitor.hpp"
65 #endif
66 #if defined(HAVE_TEUCHOSCORE_CXX11)
67 # include <type_traits>
68 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
69 #include "Teuchos_TypeTraits.hpp"
70 
71 namespace Belos {
72 
74 
75 
83  PCPGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
84  {}};
85 
91  class PCPGSolMgrOrthoFailure : public BelosError {public:
92  PCPGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
93  {}};
94 
101  class PCPGSolMgrLAPACKFailure : public BelosError {public:
102  PCPGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
103  {}};
104 
111  class PCPGSolMgrRecyclingFailure : public BelosError {public:
112  PCPGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg)
113  {}};
114 
116 
117 
141 
142  // Partial specialization for complex ScalarType.
143  // This contains a trivial implementation.
144  // See discussion in the class documentation above.
145  //
146  // FIXME (mfh 09 Sep 2015) This also is a stub for types other than
147  // float or double.
148  template<class ScalarType, class MV, class OP,
149  const bool supportsScalarType =
151  ! Teuchos::ScalarTraits<ScalarType>::isComplex>
152  class PCPGSolMgr :
153  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
154  Belos::Details::LapackSupportsScalar<ScalarType>::value &&
155  ! Teuchos::ScalarTraits<ScalarType>::isComplex>
156  {
157  static const bool scalarTypeIsSupported =
159  ! Teuchos::ScalarTraits<ScalarType>::isComplex;
160  typedef Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
162 
163  public:
165  base_type ()
166  {}
167  PCPGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
168  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
169  base_type ()
170  {}
171  virtual ~PCPGSolMgr () {}
172 
174  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
175  return Teuchos::rcp(new PCPGSolMgr<ScalarType,MV,OP,supportsScalarType>);
176  }
177  };
178 
179  template<class ScalarType, class MV, class OP>
180  class PCPGSolMgr<ScalarType, MV, OP, true> :
181  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, true> {
182  private:
185  typedef Teuchos::ScalarTraits<ScalarType> SCT;
186  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
187  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
188 
189  public:
191 
192 
199  PCPGSolMgr();
200 
236  PCPGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
237  const Teuchos::RCP<Teuchos::ParameterList> &pl );
238 
240  virtual ~PCPGSolMgr() {};
241 
243  virtual Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const {
244  return Teuchos::rcp(new PCPGSolMgr<ScalarType,MV,OP>);
245  }
247 
249 
250 
254  return *problem_;
255  }
256 
259  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
260 
263  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
264 
270  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
271  return Teuchos::tuple(timerSolve_);
272  }
273 
280  return achievedTol_;
281  }
282 
284  int getNumIters() const {
285  return numIters_;
286  }
287 
290  bool isLOADetected() const { return false; }
291 
293 
295 
296 
298  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
299 
301  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
302 
304 
306 
307 
311  void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
313 
315 
316 
334  ReturnType solve();
335 
337 
340 
342  std::string description() const;
343 
345 
346  private:
347 
348  // In the A-inner product, perform an RRQR decomposition without using A unless absolutely necessary. Given
349  // the seed space U and C = A U, find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU.
350  int ARRQR(int numVecs, int numOrthVecs, const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
351 
352  // Linear problem.
353  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
354 
355  // Output manager.
356  Teuchos::RCP<OutputManager<ScalarType> > printer_;
357  Teuchos::RCP<std::ostream> outputStream_;
358 
359  // Status test.
360  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
361  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
362  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
363  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
364 
365  // Orthogonalization manager.
366  Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
367 
368  // Current parameter list.
369  Teuchos::RCP<Teuchos::ParameterList> params_;
370 
371  // Default solver values.
372  static constexpr int maxIters_default_ = 1000;
373  static constexpr int deflatedBlocks_default_ = 2;
374  static constexpr int savedBlocks_default_ = 16;
375  static constexpr int verbosity_default_ = Belos::Errors;
376  static constexpr int outputStyle_default_ = Belos::General;
377  static constexpr int outputFreq_default_ = -1;
378  static constexpr const char * label_default_ = "Belos";
379  static constexpr const char * orthoType_default_ = "ICGS";
380 // https://stackoverflow.com/questions/24398102/constexpr-and-initialization-of-a-static-const-void-pointer-with-reinterpret-cas
381 #if defined(_WIN32) && defined(__clang__)
382  static constexpr std::ostream * outputStream_default_ =
383  __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
384 #else
385  static constexpr std::ostream * outputStream_default_ = &std::cout;
386 #endif
387 
388  //
389  // Current solver values.
390  //
391 
394 
397 
400 
403 
406 
407  int deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
408  std::string orthoType_;
409 
410  // Recycled subspace, its image and the residual
411  Teuchos::RCP<MV> U_, C_, R_;
412 
413  // Actual dimension of current recycling subspace (<= savedBlocks_ )
414  int dimU_;
415 
416  // Timers.
417  std::string label_;
418  Teuchos::RCP<Teuchos::Time> timerSolve_;
419 
420  // Internal state variables.
421  bool isSet_;
422  };
423 
424 
425 // Empty Constructor
426 template<class ScalarType, class MV, class OP>
428  outputStream_(Teuchos::rcp(outputStream_default_,false)),
429  convtol_(DefaultSolverParameters::convTol),
430  orthoKappa_(DefaultSolverParameters::orthoKappa),
431  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
432  numIters_(0),
433  maxIters_(maxIters_default_),
434  deflatedBlocks_(deflatedBlocks_default_),
435  savedBlocks_(savedBlocks_default_),
436  verbosity_(verbosity_default_),
437  outputStyle_(outputStyle_default_),
438  outputFreq_(outputFreq_default_),
439  orthoType_(orthoType_default_),
440  dimU_(0),
441  label_(label_default_),
442  isSet_(false)
443 {}
444 
445 
446 // Basic Constructor
447 template<class ScalarType, class MV, class OP>
449  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
450  const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
451  problem_(problem),
452  outputStream_(Teuchos::rcp(outputStream_default_,false)),
453 
454  convtol_(DefaultSolverParameters::convTol),
455  orthoKappa_(DefaultSolverParameters::orthoKappa),
456  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
457  numIters_(0),
458  maxIters_(maxIters_default_),
459  deflatedBlocks_(deflatedBlocks_default_),
460  savedBlocks_(savedBlocks_default_),
461  verbosity_(verbosity_default_),
462  outputStyle_(outputStyle_default_),
463  outputFreq_(outputFreq_default_),
464  orthoType_(orthoType_default_),
465  dimU_(0),
466  label_(label_default_),
467  isSet_(false)
468 {
469  TEUCHOS_TEST_FOR_EXCEPTION(
470  problem_.is_null (), std::invalid_argument,
471  "Belos::PCPGSolMgr two-argument constructor: "
472  "'problem' is null. You must supply a non-null Belos::LinearProblem "
473  "instance when calling this constructor.");
474 
475  if (! pl.is_null ()) {
476  // Set the parameters using the list that was passed in.
477  setParameters (pl);
478  }
479 }
480 
481 
482 template<class ScalarType, class MV, class OP>
483 void PCPGSolMgr<ScalarType,MV,OP,true>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
484 {
485  // Create the internal parameter list if ones doesn't already exist.
486  if (params_ == Teuchos::null) {
487  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
488  }
489  else {
490  params->validateParameters(*getValidParameters());
491  }
492 
493  // Check for maximum number of iterations
494  if (params->isParameter("Maximum Iterations")) {
495  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
496 
497  // Update parameter in our list and in status test.
498  params_->set("Maximum Iterations", maxIters_);
499  if (maxIterTest_!=Teuchos::null)
500  maxIterTest_->setMaxIters( maxIters_ );
501  }
502 
503  // Check for the maximum numbers of saved and deflated blocks.
504  if (params->isParameter("Num Saved Blocks")) {
505  savedBlocks_ = params->get("Num Saved Blocks",savedBlocks_default_);
506  TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
507  "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
508 
509  // savedBlocks > number of matrix rows and columns, not known in parameters.
510  //TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ >= maxIters_, std::invalid_argument,
511  //"Belos::PCPGSolMgr: \"Num Saved Blocks\" must be less than \"Maximum Iterations\".");
512 
513  // Update parameter in our list.
514  params_->set("Num Saved Blocks", savedBlocks_);
515  }
516  if (params->isParameter("Num Deflated Blocks")) {
517  deflatedBlocks_ = params->get("Num Deflated Blocks",deflatedBlocks_default_);
518  TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
519  "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
520 
521  TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
522  "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
523 
524  // Update parameter in our list.
525  // The static_cast is for clang link issues with the constexpr before c++17
526  params_->set("Num Deflated Blocks", static_cast<int>(deflatedBlocks_));
527  }
528 
529  // Check to see if the timer label changed.
530  if (params->isParameter("Timer Label")) {
531  std::string tempLabel = params->get("Timer Label", label_default_);
532 
533  // Update parameter in our list and solver timer
534  if (tempLabel != label_) {
535  label_ = tempLabel;
536  params_->set("Timer Label", label_);
537  std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
538 #ifdef BELOS_TEUCHOS_TIME_MONITOR
539  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
540 #endif
541  if (ortho_ != Teuchos::null) {
542  ortho_->setLabel( label_ );
543  }
544  }
545  }
546 
547  // Check for a change in verbosity level
548  if (params->isParameter("Verbosity")) {
549  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
550  verbosity_ = params->get("Verbosity", verbosity_default_);
551  } else {
552  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
553  }
554 
555  // Update parameter in our list.
556  params_->set("Verbosity", verbosity_);
557  if (printer_ != Teuchos::null)
558  printer_->setVerbosity(verbosity_);
559  }
560 
561  // Check for a change in output style
562  if (params->isParameter("Output Style")) {
563  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
564  outputStyle_ = params->get("Output Style", outputStyle_default_);
565  } else {
566  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
567  }
568 
569  // Reconstruct the convergence test if the explicit residual test is not being used.
570  params_->set("Output Style", outputStyle_);
571  outputTest_ = Teuchos::null;
572  }
573 
574  // output stream
575  if (params->isParameter("Output Stream")) {
576  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
577 
578  // Update parameter in our list.
579  params_->set("Output Stream", outputStream_);
580  if (printer_ != Teuchos::null)
581  printer_->setOStream( outputStream_ );
582  }
583 
584  // frequency level
585  if (verbosity_ & Belos::StatusTestDetails) {
586  if (params->isParameter("Output Frequency")) {
587  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
588  }
589 
590  // Update parameter in out list and output status test.
591  params_->set("Output Frequency", outputFreq_);
592  if (outputTest_ != Teuchos::null)
593  outputTest_->setOutputFrequency( outputFreq_ );
594  }
595 
596  // Create output manager if we need to.
597  if (printer_ == Teuchos::null) {
598  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
599  }
600 
601  // Check if the orthogonalization changed.
602  bool changedOrthoType = false;
603  if (params->isParameter("Orthogonalization")) {
604  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
605  if (tempOrthoType != orthoType_) {
606  orthoType_ = tempOrthoType;
607  changedOrthoType = true;
608  }
609  }
610  params_->set("Orthogonalization", orthoType_);
611 
612  // Check which orthogonalization constant to use.
613  if (params->isParameter("Orthogonalization Constant")) {
614  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
615  orthoKappa_ = params->get ("Orthogonalization Constant",
616  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
617  }
618  else {
619  orthoKappa_ = params->get ("Orthogonalization Constant",
621  }
622 
623  // Update parameter in our list.
624  params_->set("Orthogonalization Constant",orthoKappa_);
625  if (orthoType_=="DGKS") {
626  if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
627  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
628  }
629  }
630  }
631 
632  // Create orthogonalization manager if we need to.
633  if (ortho_ == Teuchos::null || changedOrthoType) {
635  Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
636  if (orthoType_=="DGKS" && orthoKappa_ > 0) {
637  paramsOrtho->set ("depTol", orthoKappa_ );
638  }
639 
640  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
641  }
642 
643  // Convergence
644  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
645  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
646 
647  // Check for convergence tolerance
648  if (params->isParameter("Convergence Tolerance")) {
649  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
650  convtol_ = params->get ("Convergence Tolerance",
651  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
652  }
653  else {
654  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
655  }
656 
657  // Update parameter in our list and residual tests.
658  params_->set("Convergence Tolerance", convtol_);
659  if (convTest_ != Teuchos::null)
660  convTest_->setTolerance( convtol_ );
661  }
662 
663  // Create status tests if we need to.
664 
665  // Basic test checks maximum iterations and native residual.
666  if (maxIterTest_ == Teuchos::null)
667  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
668 
669  if (convTest_ == Teuchos::null)
670  convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
671 
672  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
673 
674  // Create the status test output class.
675  // This class manages and formats the output from the status test.
676  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
677  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
678 
679  // Set the solver string for the output test
680  std::string solverDesc = " PCPG ";
681  outputTest_->setSolverDesc( solverDesc );
682 
683  // Create the timer if we need to.
684  if (timerSolve_ == Teuchos::null) {
685  std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
686 #ifdef BELOS_TEUCHOS_TIME_MONITOR
687  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
688 #endif
689  }
690 
691  // Inform the solver manager that the current parameters were set.
692  isSet_ = true;
693 }
694 
695 
696 template<class ScalarType, class MV, class OP>
697 Teuchos::RCP<const Teuchos::ParameterList>
699 {
700  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
701  if (is_null(validPL)) {
702  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
703  // Set all the valid parameters and their default values.
704  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
705  "The relative residual tolerance that needs to be achieved by the\n"
706  "iterative solver in order for the linear system to be declared converged.");
707  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
708  "The maximum number of iterations allowed for each\n"
709  "set of RHS solved.");
710  pl->set("Num Deflated Blocks", static_cast<int>(deflatedBlocks_default_),
711  "The maximum number of vectors in the seed subspace." );
712  pl->set("Num Saved Blocks", static_cast<int>(savedBlocks_default_),
713  "The maximum number of vectors saved from old Krylov subspaces." );
714  pl->set("Verbosity", static_cast<int>(verbosity_default_),
715  "What type(s) of solver information should be outputted\n"
716  "to the output stream.");
717  pl->set("Output Style", static_cast<int>(outputStyle_default_),
718  "What style is used for the solver information outputted\n"
719  "to the output stream.");
720  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
721  "How often convergence information should be outputted\n"
722  "to the output stream.");
723  pl->set("Output Stream", Teuchos::rcp(outputStream_default_,false),
724  "A reference-counted pointer to the output stream where all\n"
725  "solver output is sent.");
726  pl->set("Timer Label", static_cast<const char *>(label_default_),
727  "The string to use as a prefix for the timer labels.");
728  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
729  "The type of orthogonalization to use: DGKS, ICGS, IMGS");
730  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
731  "The constant used by DGKS orthogonalization to determine\n"
732  "whether another step of classical Gram-Schmidt is necessary.");
733  validPL = pl;
734  }
735  return validPL;
736 }
737 
738 
739 // solve()
740 template<class ScalarType, class MV, class OP>
742 
743  // Set the current parameters if are not set already.
744  if (!isSet_) { setParameters( params_ ); }
745 
746  Teuchos::LAPACK<int,ScalarType> lapack;
747  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
748  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
749 
750  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,PCPGSolMgrLinearProblemFailure,
751  "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
752 
753  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PCPGSolMgrLinearProblemFailure,
754  "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
755 
756  // Create indices for the linear systems to be solved.
757  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
758  std::vector<int> currIdx(1);
759  currIdx[0] = 0;
760 
761  bool debug = false;
762 
763  // Inform the linear problem of the current linear system to solve.
764  problem_->setLSIndex( currIdx ); // block size == 1
765 
766  // Assume convergence is achieved, then let any failed convergence set this to false.
767  bool isConverged = true;
768 
770  // PCPG iteration parameter list
771  Teuchos::ParameterList plist;
772  plist.set("Saved Blocks", savedBlocks_);
773  plist.set("Block Size", 1);
774  plist.set("Keep Diagonal", true);
775  plist.set("Initialize Diagonal", true);
776 
778  // PCPG solver
779 
780  Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
781  pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
782  // Number of iterations required to generate initial recycle space (if needed)
783 
784  // Enter solve() iterations
785  {
786 #ifdef BELOS_TEUCHOS_TIME_MONITOR
787  Teuchos::TimeMonitor slvtimer(*timerSolve_);
788 #endif
789  while ( numRHS2Solve > 0 ) { // test for quick return
790 
791  // Reset the status test.
792  outputTest_->reset();
793 
794  // Create the first block in the current Krylov basis (residual).
795  if (R_ == Teuchos::null)
796  R_ = MVT::Clone( *(problem_->getRHS()), 1 );
797 
798  problem_->computeCurrResVec( &*R_ );
799 
800 
801  // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
802  // TODO: ensure hypothesis right here ... I have to think about use cases.
803 
804  if( U_ != Teuchos::null ){
805  // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
806 
807  // possibly over solved equation ... I want residual norms
808  // relative to the initial residual, not what I am about to compute.
809  Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
810  std::vector<MagnitudeType> rnorm0(1);
811  MVT::MvNorm( *R_, rnorm0 ); // rnorm0 = norm(R_);
812 
813  // Z := U_'*R_; xo += U_*Z ;R_ -= C_*Z
814  std::cout << "Solver Manager: dimU_ = " << dimU_ << std::endl;
815  Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
816 
817  Teuchos::RCP<const MV> Uactive, Cactive;
818  std::vector<int> active_columns( dimU_ );
819  for (int i=0; i < dimU_; ++i) active_columns[i] = i;
820  Uactive = MVT::CloneView(*U_, active_columns);
821  Cactive = MVT::CloneView(*C_, active_columns);
822 
823  if( debug ){
824  std::cout << " Solver Manager : check duality of seed basis " << std::endl;
825  Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
826  MVT::MvTransMv( one, *Uactive, *Cactive, H );
827  H.print( std::cout );
828  }
829 
830  MVT::MvTransMv( one, *Uactive, *R_, Z );
831  Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
832  MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU ); // UZ
833  MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec ); // xo += tmp;
834  MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU ); // CZ
835  MVT::MvAddMv( -one, *tempU, one, *R_, *R_ ); // R_ -= tmp;
836  std::vector<MagnitudeType> rnorm(1);
837  MVT::MvNorm( *R_, rnorm );
838  if( rnorm[0] < rnorm0[0] * .001 ){ //reorthogonalize
839  MVT::MvTransMv( one, *Uactive, *R_, Z );
840  MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
841  MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec ); // xo += UZ;
842  MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
843  MVT::MvAddMv( -one, *tempU, one, *R_, *R_ ); // R_ -= CZ;
844  }
845  Uactive = Teuchos::null;
846  Cactive = Teuchos::null;
847  tempU = Teuchos::null;
848  }
849  else {
850  dimU_ = 0;
851  }
852 
853 
854  // Set the new state and initialize the solver.
855  PCPGIterState<ScalarType,MV> pcpgState; // fails if R == null.
856 
857  pcpgState.R = R_;
858  if( U_ != Teuchos::null ) pcpgState.U = U_;
859  if( C_ != Teuchos::null ) pcpgState.C = C_;
860  if( dimU_ > 0 ) pcpgState.curDim = dimU_;
861  pcpg_iter->initialize(pcpgState);
862 
863  // treat initialize() exceptions here? how to use try-catch-throw? DMD
864 
865  // Get the current number of deflated blocks with the PCPG iteration
866  dimU_ = pcpgState.curDim;
867  if( !dimU_ ) printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
868  pcpg_iter->resetNumIters();
869 
870  if( dimU_ > savedBlocks_ )
871  std::cout << "Error: dimU_ = " << dimU_ << " > savedBlocks_ = " << savedBlocks_ << std::endl;
872 
873  while(1) { // dummy loop for break
874 
875  // tell pcpg_iter to iterate
876  try {
877  if( debug ) printf("********** Calling iterate...\n");
878  pcpg_iter->iterate();
879 
881  //
882  // check convergence first
883  //
885  if ( convTest_->getStatus() == Passed ) {
886  // we have convergence
887  break; // break from while(1){pcpg_iter->iterate()}
888  }
890  //
891  // check for maximum iterations
892  //
894  else if ( maxIterTest_->getStatus() == Passed ) {
895  // we don't have convergence
896  isConverged = false;
897  break; // break from while(1){pcpg_iter->iterate()}
898  }
899  else {
900 
902  //
903  // we returned from iterate(), but none of our status tests Passed.
904  // Something is wrong, and it is probably the developers fault.
905  //
907 
908  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
909  "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
910  } // end if
911  } // end try
912  catch (const PCPGIterOrthoFailure &e) {
913 
914  // Check to see if the most recent solution yielded convergence.
915  sTest_->checkStatus( &*pcpg_iter );
916  if (convTest_->getStatus() != Passed)
917  isConverged = false;
918  break;
919  }
920  catch (const std::exception &e) {
921  printer_->stream(Errors) << "Error! Caught exception in PCPGIter::iterate() at iteration "
922  << pcpg_iter->getNumIters() << std::endl
923  << e.what() << std::endl;
924  throw;
925  }
926  } // end of while(1)
927 
928  // Update the linear problem.
929  Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
930  problem_->updateSolution( update, true );
931 
932  // Inform the linear problem that we are finished with this block linear system.
933  problem_->setCurrLS();
934 
935  // Get the state. How did pcpgState die?
936  PCPGIterState<ScalarType,MV> oldState = pcpg_iter->getState();
937 
938  dimU_ = oldState.curDim;
939  int q = oldState.prevUdim;
940 
941  std::cout << "SolverManager: dimU_ " << dimU_ << " prevUdim= " << q << std::endl;
942 
943  if( q > deflatedBlocks_ )
944  std::cout << "SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
945 
946  int rank;
947  if( dimU_ > q ){ // Orthogonalize [U;C](:,prevUdim:dimU_)
948  //Given the seed space U and C = A U for some symmetric positive definite A,
949  //find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU
950 
951  //oldState.D->print( std::cout ); D = diag( C'*U )
952 
953  U_ = oldState.U; //MVT::MvPrint( *U_, std::cout );
954  C_ = oldState.C; //MVT::MvPrint( *C_, std::cout );
955  rank = ARRQR(dimU_,q, *oldState.D );
956  if( rank < dimU_ ) {
957  std::cout << " rank decreased in ARRQR, something to do? " << std::endl;
958  }
959  dimU_ = rank;
960 
961  } // Now U_ and C_ = AU are dual bases.
962 
963  if( dimU_ > deflatedBlocks_ ){
964 
965  if( !deflatedBlocks_ ){
966  U_ = Teuchos::null;
967  C_ = Teuchos::null;
968  dimU_ = deflatedBlocks_;
969  break;
970  }
971 
972  bool Harmonic = false; // (Harmonic) Ritz vectors
973 
974  Teuchos::RCP<MV> Uorth;
975 
976  std::vector<int> active_cols( dimU_ );
977  for (int i=0; i < dimU_; ++i) active_cols[i] = i;
978 
979  if( Harmonic ){
980  Uorth = MVT::CloneCopy(*C_, active_cols);
981  }
982  else{
983  Uorth = MVT::CloneCopy(*U_, active_cols);
984  }
985 
986  // Explicitly construct Q and R factors
987  Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
988  rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,false));
989  Uorth = Teuchos::null;
990  // TODO: During the previous solve, the matrix that normalizes U(1:q) was computed and discarded.
991  // One might save it, reuse it here, and just normalize columns U(q+1:dimU_) here.
992 
993  // throw an error if U is both A-orthonormal and rank deficient
994  TEUCHOS_TEST_FOR_EXCEPTION(rank < dimU_,PCPGSolMgrOrthoFailure,
995  "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
996 
997 
998  // R VT' = Ur S,
999  Teuchos::SerialDenseMatrix<int,ScalarType> VT; // Not referenced
1000  Teuchos::SerialDenseMatrix<int,ScalarType> Ur; // Not referenced
1001  int lwork = 5*dimU_; // minimal, extra computation < 67*dimU_
1002  int info = 0; // Hermite
1003  int lrwork = 1;
1004  if( problem_->isHermitian() ) lrwork = dimU_;
1005  std::vector<ScalarType> work(lwork); //
1006  std::vector<ScalarType> Svec(dimU_); //
1007  std::vector<ScalarType> rwork(lrwork);
1008  lapack.GESVD('N', 'O',
1009  R.numRows(),R.numCols(),R.values(), R.numRows(),
1010  &Svec[0],
1011  Ur.values(),1,
1012  VT.values(),1, // Output: VT stored in R
1013  &work[0], lwork,
1014  &rwork[0], &info);
1015 
1016  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, PCPGSolMgrLAPACKFailure,
1017  "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
1018 
1019  if( work[0] != 67. * dimU_ )
1020  std::cout << " SVD " << dimU_ << " lwork " << work[0] << std::endl;
1021  for( int i=0; i< dimU_; i++)
1022  std::cout << i << " " << Svec[i] << std::endl;
1023 
1024  Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R, Teuchos::TRANS);
1025 
1026  int startRow = 0, startCol = 0;
1027  if( Harmonic )
1028  startCol = dimU_ - deflatedBlocks_;
1029 
1030  Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
1031  wholeV,
1032  wholeV.numRows(),
1033  deflatedBlocks_,
1034  startRow,
1035  startCol);
1036  std::vector<int> active_columns( dimU_ );
1037  std::vector<int> def_cols( deflatedBlocks_ );
1038  for (int i=0; i < dimU_; ++i) active_columns[i] = i;
1039  for (int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
1040 
1041  Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
1042  Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
1043  MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive ); // U:= U*V
1044  Ucopy = Teuchos::null;
1045  Uactive = Teuchos::null;
1046  Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
1047  Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
1048  MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive ); // C:= C*V
1049  Ccopy = Teuchos::null;
1050  Cactive = Teuchos::null;
1051  dimU_ = deflatedBlocks_;
1052  }
1053  printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << dimU_ << std::endl << std::endl;
1054 
1055  // Inform the linear problem that we are finished with this block linear system.
1056  problem_->setCurrLS();
1057 
1058  // Update indices for the linear systems to be solved.
1059  numRHS2Solve -= 1;
1060  if ( numRHS2Solve > 0 ) {
1061  currIdx[0]++;
1062 
1063  // Set the next indices.
1064  problem_->setLSIndex( currIdx );
1065  }
1066  else {
1067  currIdx.resize( numRHS2Solve );
1068  }
1069  }// while ( numRHS2Solve > 0 )
1070  }
1071 
1072  // print final summary
1073  sTest_->print( printer_->stream(FinalSummary) );
1074 
1075  // print timing information
1076 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1077  // Calling summarize() can be expensive, so don't call unless the
1078  // user wants to print out timing details. summarize() will do all
1079  // the work even if it's passed a "black hole" output stream.
1080  if (verbosity_ & TimingDetails)
1081  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1082 #endif
1083 
1084  // Save the convergence test value ("achieved tolerance") for this solve.
1085  {
1086  using Teuchos::rcp_dynamic_cast;
1087  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1088  // testValues is nonnull and not persistent.
1089  const std::vector<MagnitudeType>* pTestValues =
1090  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1091 
1092  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1093  "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1094  "method returned NULL. Please report this bug to the Belos developers.");
1095 
1096  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1097  "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1098  "method returned a vector of length zero. Please report this bug to the "
1099  "Belos developers.");
1100 
1101  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1102  // achieved tolerances for all vectors in the current solve(), or
1103  // just for the vectors from the last deflation?
1104  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1105  }
1106 
1107  // get iteration information for this solve
1108  numIters_ = maxIterTest_->getNumIters();
1109 
1110  if (!isConverged) {
1111  return Unconverged; // return from PCPGSolMgr::solve()
1112  }
1113  return Converged; // return from PCPGSolMgr::solve()
1114 }
1115 
1116 // A-orthogonalize the Seed Space
1117 // Note that Anasazi::GenOrthoManager provides simplified versions of the algorithm,
1118 // that are not rank revealing, and are not designed for PCPG in other ways too.
1119 template<class ScalarType, class MV, class OP>
1120 int PCPGSolMgr<ScalarType,MV,OP,true>::ARRQR(int p, int q, const Teuchos::SerialDenseMatrix<int,ScalarType>& D)
1121 {
1122  using Teuchos::RCP;
1123  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1124  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1125 
1126  // Allocate memory for scalars.
1127  Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
1128  Teuchos::SerialDenseMatrix<int,ScalarType> gamma( 1, 1 );
1129  Teuchos::SerialDenseMatrix<int,ScalarType> anorm( 1, 1 );
1130  std::vector<int> curind(1);
1131  std::vector<int> ipiv(p - q); // RRQR Pivot indices
1132  std::vector<ScalarType> Pivots(p); // RRQR Pivots
1133  int i, imax, j, l;
1134  ScalarType rteps = 1.5e-8;
1135 
1136  // Scale such that diag( U'C) = I
1137  for( i = q ; i < p ; i++ ){
1138  ipiv[i-q] = i;
1139  curind[0] = i;
1140  RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1141  RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1142  anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
1143  MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1144  MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1145  Pivots[i] = one;
1146  }
1147 
1148  for( i = q ; i < p ; i++ ){
1149  if( q < i && i < p-1 ){ // Find the largest pivot
1150  imax = i;
1151  l = ipiv[imax-q];
1152  for( j = i+1 ; j < p ; j++ ){
1153  const int k = ipiv[j-q];
1154  if( Pivots[k] > Pivots[l] ){
1155  imax = j;
1156  l = k;
1157  }
1158  } // end for
1159  if( imax > i ){
1160  l = ipiv[imax-q]; // swap ipiv( imax ) and ipiv(i+1)
1161  ipiv[imax-q] = ipiv[i-q];
1162  ipiv[i-q] = l;
1163  }
1164  } // largest pivot found
1165  int k = ipiv[i-q];
1166 
1167  if( Pivots[k] > 1.5625e-2 ){
1168  anorm(0,0) = Pivots[k]; // A-norm of u
1169  }
1170  else{ // anorm(0,0) = sqrt( U(:,k)'*C(:,k) );
1171  curind[0] = k;
1172  RCP<const MV> P = MVT::CloneView(*U_,curind);
1173  RCP<const MV> AP = MVT::CloneView(*C_,curind);
1174  MVT::MvTransMv( one, *P, *AP, anorm );
1175  anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
1176  }
1177  if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
1178  /*
1179  C(:,k) = A*U(:,k); % Change C
1180  fixC = U(:, ipiv(1:i-1) )'*C(:,k);
1181  U(:,k) = U(:,k) - U(:, ipiv(1:i-1) )*fixC;
1182  C(:,k) = C(:,k) - C(:, ipiv(1:i-1) )*fixC;
1183  anorm = sqrt( U(:,k)'*C(:,k) );
1184  */
1185  std::cout << "ARRQR: Bad case not implemented" << std::endl;
1186  }
1187  if( anorm(0,0) < rteps ){ // rank [U;C] = i-1
1188  std::cout << "ARRQR : deficient case not implemented " << std::endl;
1189  //U = U(:, ipiv(1:i-1) );
1190  //C = C(:, ipiv(1:i-1) );
1191  p = q + i;
1192  // update curDim_ in State
1193  break;
1194  }
1195  curind[0] = k;
1196  RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1197  RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1198  MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P ); // U(:,k) = U(:,k)/anorm;
1199  MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP ); // C(:,k) = C(:,k)/anorm;
1200  P = Teuchos::null;
1201  AP = Teuchos::null;
1202  Pivots[k] = one; // delete, for diagonostic purposes
1203  P = MVT::CloneViewNonConst(*U_,curind); // U(:,k)
1204  AP = MVT::CloneViewNonConst(*C_,curind); // C(:,k)
1205  for( j = i+1 ; j < p ; j++ ){
1206  l = ipiv[j-q]; // ahhh
1207  curind[0] = l;
1208  RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind); // segmentation fault, j=i+1=5
1209  MVT::MvTransMv( one, *Q, *AP, alpha); // alpha(0,0) = U(:,l)'*C(:,k);
1210  MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q ); // U(:,l) -= U(:,k) * alpha(0,0);
1211  Q = Teuchos::null;
1212  RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
1213  MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ ); // C(:,l) -= C(:,l) - C(:,k) * alpha(0,0);
1214  AQ = Teuchos::null;
1215  gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
1216  if( gamma(0,0) > 0){
1217  Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
1218  }
1219  else {
1220  Pivots[l] = zero; //rank deficiency revealed
1221  }
1222  }
1223  }
1224  return p;
1225 }
1226 
1227 // The method returns a string describing the solver manager.
1228 template<class ScalarType, class MV, class OP>
1230 {
1231  std::ostringstream oss;
1232  oss << "Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1233  oss << "{";
1234  oss << "Ortho Type='"<<orthoType_;
1235  oss << "}";
1236  return oss.str();
1237 }
1238 
1239 } // end Belos namespace
1240 
1241 #endif /* BELOS_PCPG_SOLMGR_HPP */
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:299
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Collection of types and exceptions used within the Belos solvers.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
PCPGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Class which manages the output and verbosity of the Belos solvers.
static const bool scalarTypeIsSupported
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Belos concrete class to iterate Preconditioned Conjugate Projected Gradients.
int numIters_
Number of iterations taken by the last solve() invocation.
PCPGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonormal...
Teuchos::RCP< MV > R
The current residual.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
PCPGSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspace...
A factory class for generating StatusTestOutput objects.
int curDim
The current dimension of the reduction.
An implementation of StatusTestResNorm using a family of residual norms.
PCPGSolMgrOrthoFailure(const std::string &what_arg)
int maxIters_
Maximum iteration count (read from parameter list).
MagnitudeType orthoKappa_
Orthogonalization parameter (read from parameter list).
Structure to contain pointers to PCPGIter state variables.
PCPGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:293
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > D
The current Hessenberg matrix.
OperatorTraits< ScalarType, MV, OP > OPT
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:206
Pure virtual base class which describes the basic interface for a solver manager. ...
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
virtual Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const
clone for Inverted Injection (DII)
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
A linear system to solve, and its associated information.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
Class which describes the linear problem to be solved by the iterative solver.
PCPGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
PCPGSolMgrLAPACKFailure(const std::string &what_arg)
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< OutputManager< ScalarType > > printer_
Teuchos::RCP< MV > U
The recycled subspace.
int setTolerance(MagnitudeType tolerance)
Set the value of the tolerance.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
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.
PCPGIterOrthoFailure is thrown when the PCPGIter object is unable to compute independent direction ve...
Teuchos::ScalarTraits< ScalarType > SCT
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
PCPG iterative linear solver.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
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.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:283
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Details::SolverManagerRequiresRealLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
PCPGSolMgrRecyclingFailure(const std::string &what_arg)
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed...
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< Teuchos::ParameterList > params_
PCPGSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
Return an instance of the specified MatOrthoManager subclass.