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

Generated on Fri Jun 5 2020 10:20:43 for Belos by doxygen 1.8.5