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