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