Belos Package Browser (Single Doxygen Collection)  Development
BelosPCPGSolMgr.hpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 #ifndef BELOS_PCPG_SOLMGR_HPP
43 #define BELOS_PCPG_SOLMGR_HPP
44 
48 
49 #include "BelosConfigDefs.hpp"
50 #include "BelosTypes.hpp"
51 #include "BelosLinearProblem.hpp"
52 #include "BelosSolverManager.hpp"
53 
54 #include "BelosPCPGIter.hpp"
55 
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,
165 
166  public:
168  base_type ()
169  {}
172  base_type ()
173  {}
174  virtual ~PCPGSolMgr () {}
175  };
176 
177  template<class ScalarType, class MV, class OP>
178  class PCPGSolMgr<ScalarType, MV, OP, true> :
179  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, true> {
180  private:
186 
187  public:
189 
190 
197  PCPGSolMgr();
198 
236 
238  virtual ~PCPGSolMgr() {};
240 
242 
243 
247  return *problem_;
248  }
249 
252  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
253 
257 
264  return Teuchos::tuple(timerSolve_);
265  }
266 
273  return achievedTol_;
274  }
275 
277  int getNumIters() const {
278  return numIters_;
279  }
280 
283  bool isLOADetected() const { return false; }
284 
286 
288 
289 
291  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
292 
294  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
295 
297 
299 
300 
304  void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
306 
308 
309 
327  ReturnType solve();
328 
330 
333 
335  std::string description() const;
336 
338 
339  private:
340 
341  // In the A-inner product, perform an RRQR decomposition without using A unless absolutely necessary. Given
342  // 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.
343  int ARRQR(int numVecs, int numOrthVecs, const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
344 
345  // Linear problem.
347 
348  // Output manager.
351 
352  // Status test.
357 
358  // Orthogonalization manager.
360 
361  // Current parameter list.
363 
364  // Default solver values.
367  static const int maxIters_default_;
368  static const int deflatedBlocks_default_;
369  static const int savedBlocks_default_;
370  static const int verbosity_default_;
371  static const int outputStyle_default_;
372  static const int outputFreq_default_;
373  static const std::string label_default_;
374  static const std::string orthoType_default_;
376 
377  //
378  // Current solver values.
379  //
380 
383 
386 
389 
392 
395 
396  int deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
397  std::string orthoType_;
398 
399  // Recycled subspace, its image and the residual
401 
402  // Actual dimension of current recycling subspace (<= savedBlocks_ )
403  int dimU_;
404 
405  // Timers.
406  std::string label_;
408 
409  // Internal state variables.
410  bool isSet_;
411  };
412 
413 
414 // Default solver values.
415 template<class ScalarType, class MV, class OP>
418 
419 template<class ScalarType, class MV, class OP>
422 
423 template<class ScalarType, class MV, class OP>
425 
426 template<class ScalarType, class MV, class OP>
428 
429 template<class ScalarType, class MV, class OP>
431 
432 template<class ScalarType, class MV, class OP>
434 
435 template<class ScalarType, class MV, class OP>
437 
438 template<class ScalarType, class MV, class OP>
440 
441 template<class ScalarType, class MV, class OP>
442 const std::string PCPGSolMgr<ScalarType,MV,OP,true>::label_default_ = "Belos";
443 
444 template<class ScalarType, class MV, class OP>
446 
447 template<class ScalarType, class MV, class OP>
449 
450 
451 // Empty Constructor
452 template<class ScalarType, class MV, class OP>
454  outputStream_(outputStream_default_),
455  convtol_(convtol_default_),
456  orthoKappa_(orthoKappa_default_),
457  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
458  numIters_(0),
459  maxIters_(maxIters_default_),
460  deflatedBlocks_(deflatedBlocks_default_),
461  savedBlocks_(savedBlocks_default_),
462  verbosity_(verbosity_default_),
463  outputStyle_(outputStyle_default_),
464  outputFreq_(outputFreq_default_),
465  orthoType_(orthoType_default_),
466  dimU_(0),
467  label_(label_default_),
468  isSet_(false)
469 {}
470 
471 
472 // Basic Constructor
473 template<class ScalarType, class MV, class OP>
477  problem_(problem),
478  outputStream_(outputStream_default_),
479  convtol_(convtol_default_),
480  orthoKappa_(orthoKappa_default_),
481  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
482  numIters_(0),
483  maxIters_(maxIters_default_),
484  deflatedBlocks_(deflatedBlocks_default_),
485  savedBlocks_(savedBlocks_default_),
486  verbosity_(verbosity_default_),
487  outputStyle_(outputStyle_default_),
488  outputFreq_(outputFreq_default_),
489  orthoType_(orthoType_default_),
490  dimU_(0),
491  label_(label_default_),
492  isSet_(false)
493 {
495  problem_.is_null (), std::invalid_argument,
496  "Belos::PCPGSolMgr two-argument constructor: "
497  "'problem' is null. You must supply a non-null Belos::LinearProblem "
498  "instance when calling this constructor.");
499 
500  if (! pl.is_null ()) {
501  // Set the parameters using the list that was passed in.
502  setParameters (pl);
503  }
504 }
505 
506 
507 template<class ScalarType, class MV, class OP>
509 {
510  // Create the internal parameter list if ones doesn't already exist.
511  if (params_ == Teuchos::null) {
512  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
513  }
514  else {
515  params->validateParameters(*getValidParameters());
516  }
517 
518  // Check for maximum number of iterations
519  if (params->isParameter("Maximum Iterations")) {
520  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
521 
522  // Update parameter in our list and in status test.
523  params_->set("Maximum Iterations", maxIters_);
524  if (maxIterTest_!=Teuchos::null)
525  maxIterTest_->setMaxIters( maxIters_ );
526  }
527 
528  // Check for the maximum numbers of saved and deflated blocks.
529  if (params->isParameter("Num Saved Blocks")) {
530  savedBlocks_ = params->get("Num Saved Blocks",savedBlocks_default_);
531  TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
532  "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
533 
534  // savedBlocks > number of matrix rows and columns, not known in parameters.
535  //TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ >= maxIters_, std::invalid_argument,
536  //"Belos::PCPGSolMgr: \"Num Saved Blocks\" must be less than \"Maximum Iterations\".");
537 
538  // Update parameter in our list.
539  params_->set("Num Saved Blocks", savedBlocks_);
540  }
541  if (params->isParameter("Num Deflated Blocks")) {
542  deflatedBlocks_ = params->get("Num Deflated Blocks",deflatedBlocks_default_);
543  TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
544  "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
545 
546  TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
547  "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
548 
549  // Update parameter in our list.
550  params_->set("Num Deflated Blocks", deflatedBlocks_);
551  }
552 
553  // Check to see if the timer label changed.
554  if (params->isParameter("Timer Label")) {
555  std::string tempLabel = params->get("Timer Label", label_default_);
556 
557  // Update parameter in our list and solver timer
558  if (tempLabel != label_) {
559  label_ = tempLabel;
560  params_->set("Timer Label", label_);
561  std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
562 #ifdef BELOS_TEUCHOS_TIME_MONITOR
563  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
564 #endif
565  if (ortho_ != Teuchos::null) {
566  ortho_->setLabel( label_ );
567  }
568  }
569  }
570 
571  // Check if the orthogonalization changed.
572  if (params->isParameter("Orthogonalization")) {
573  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
574  TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS",
575  std::invalid_argument,
576  "Belos::PCPGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
577  if (tempOrthoType != orthoType_) {
578  orthoType_ = tempOrthoType;
579  // Create orthogonalization manager
580  if (orthoType_=="DGKS") {
581  if (orthoKappa_ <= 0) {
582  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
583  }
584  else {
585  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
586  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
587  }
588  }
589  else if (orthoType_=="ICGS") {
590  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
591  }
592  else if (orthoType_=="IMGS") {
593  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
594  }
595  }
596  }
597 
598  // Check which orthogonalization constant to use.
599  if (params->isParameter("Orthogonalization Constant")) {
600  orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
601 
602  // Update parameter in our list.
603  params_->set("Orthogonalization Constant",orthoKappa_);
604  if (orthoType_=="DGKS") {
605  if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
606  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
607  }
608  }
609  }
610 
611  // Check for a change in verbosity level
612  if (params->isParameter("Verbosity")) {
613  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
614  verbosity_ = params->get("Verbosity", verbosity_default_);
615  } else {
616  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
617  }
618 
619  // Update parameter in our list.
620  params_->set("Verbosity", verbosity_);
621  if (printer_ != Teuchos::null)
622  printer_->setVerbosity(verbosity_);
623  }
624 
625  // Check for a change in output style
626  if (params->isParameter("Output Style")) {
627  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
628  outputStyle_ = params->get("Output Style", outputStyle_default_);
629  } else {
630  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
631  }
632 
633  // Reconstruct the convergence test if the explicit residual test is not being used.
634  params_->set("Output Style", outputStyle_);
635  outputTest_ = Teuchos::null;
636  }
637 
638  // output stream
639  if (params->isParameter("Output Stream")) {
640  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
641 
642  // Update parameter in our list.
643  params_->set("Output Stream", outputStream_);
644  if (printer_ != Teuchos::null)
645  printer_->setOStream( outputStream_ );
646  }
647 
648  // frequency level
649  if (verbosity_ & Belos::StatusTestDetails) {
650  if (params->isParameter("Output Frequency")) {
651  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
652  }
653 
654  // Update parameter in out list and output status test.
655  params_->set("Output Frequency", outputFreq_);
656  if (outputTest_ != Teuchos::null)
657  outputTest_->setOutputFrequency( outputFreq_ );
658  }
659 
660  // Create output manager if we need to.
661  if (printer_ == Teuchos::null) {
662  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
663  }
664 
665  // Convergence
666  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
667  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
668 
669  // Check for convergence tolerance
670  if (params->isParameter("Convergence Tolerance")) {
671  convtol_ = params->get("Convergence Tolerance",convtol_default_);
672 
673  // Update parameter in our list and residual tests.
674  params_->set("Convergence Tolerance", convtol_);
675  if (convTest_ != Teuchos::null)
676  convTest_->setTolerance( convtol_ );
677  }
678 
679  // Create status tests if we need to.
680 
681  // Basic test checks maximum iterations and native residual.
682  if (maxIterTest_ == Teuchos::null)
683  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
684 
685  if (convTest_ == Teuchos::null)
686  convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
687 
688  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
689 
690  // Create the status test output class.
691  // This class manages and formats the output from the status test.
692  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
693  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
694 
695  // Set the solver string for the output test
696  std::string solverDesc = " PCPG ";
697  outputTest_->setSolverDesc( solverDesc );
698 
699 
700  // Create orthogonalization manager if we need to.
701  if (ortho_ == Teuchos::null) {
702  if (orthoType_=="DGKS") {
703  if (orthoKappa_ <= 0) {
704  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
705  }
706  else {
707  ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
708  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
709  }
710  }
711  else if (orthoType_=="ICGS") {
712  ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
713  }
714  else if (orthoType_=="IMGS") {
715  ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
716  }
717  else {
718  TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
719  "Belos::PCPGSolMgr(): Invalid orthogonalization type.");
720  }
721  }
722 
723  // Create the timer if we need to.
724  if (timerSolve_ == Teuchos::null) {
725  std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
726 #ifdef BELOS_TEUCHOS_TIME_MONITOR
727  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
728 #endif
729  }
730 
731  // Inform the solver manager that the current parameters were set.
732  isSet_ = true;
733 }
734 
735 
736 template<class ScalarType, class MV, class OP>
739 {
741  if (is_null(validPL)) {
742  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
743  // Set all the valid parameters and their default values.
744  pl->set("Convergence Tolerance", convtol_default_,
745  "The relative residual tolerance that needs to be achieved by the\n"
746  "iterative solver in order for the linear system to be declared converged.");
747  pl->set("Maximum Iterations", maxIters_default_,
748  "The maximum number of iterations allowed for each\n"
749  "set of RHS solved.");
750  pl->set("Num Deflated Blocks", deflatedBlocks_default_,
751  "The maximum number of vectors in the seed subspace." );
752  pl->set("Num Saved Blocks", savedBlocks_default_,
753  "The maximum number of vectors saved from old Krylov subspaces." );
754  pl->set("Verbosity", verbosity_default_,
755  "What type(s) of solver information should be outputted\n"
756  "to the output stream.");
757  pl->set("Output Style", outputStyle_default_,
758  "What style is used for the solver information outputted\n"
759  "to the output stream.");
760  pl->set("Output Frequency", outputFreq_default_,
761  "How often convergence information should be outputted\n"
762  "to the output stream.");
763  pl->set("Output Stream", outputStream_default_,
764  "A reference-counted pointer to the output stream where all\n"
765  "solver output is sent.");
766  pl->set("Timer Label", label_default_,
767  "The string to use as a prefix for the timer labels.");
768  // pl->set("Restart Timers", restartTimers_);
769  pl->set("Orthogonalization", orthoType_default_,
770  "The type of orthogonalization to use: DGKS, ICGS, IMGS");
771  pl->set("Orthogonalization Constant",orthoKappa_default_,
772  "The constant used by DGKS orthogonalization to determine\n"
773  "whether another step of classical Gram-Schmidt is necessary.");
774  validPL = pl;
775  }
776  return validPL;
777 }
778 
779 
780 // solve()
781 template<class ScalarType, class MV, class OP>
783 
784  // Set the current parameters if are not set already.
785  if (!isSet_) { setParameters( params_ ); }
786 
789  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
790  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
791 
793  "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
794 
796  "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
797 
798  // Create indices for the linear systems to be solved.
799  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
800  std::vector<int> currIdx(1);
801  currIdx[0] = 0;
802 
803  bool debug = false;
804 
805  // Inform the linear problem of the current linear system to solve.
806  problem_->setLSIndex( currIdx ); // block size == 1
807 
808  // Assume convergence is achieved, then let any failed convergence set this to false.
809  bool isConverged = true;
810 
812  // PCPG iteration parameter list
814  plist.set("Saved Blocks", savedBlocks_);
815  plist.set("Block Size", 1);
816  plist.set("Keep Diagonal", true);
817  plist.set("Initialize Diagonal", true);
818 
820  // PCPG solver
821 
823  pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
824  // Number of iterations required to generate initial recycle space (if needed)
825 
826  // Enter solve() iterations
827  {
828 #ifdef BELOS_TEUCHOS_TIME_MONITOR
829  Teuchos::TimeMonitor slvtimer(*timerSolve_);
830 #endif
831  while ( numRHS2Solve > 0 ) { // test for quick return
832 
833  // Reset the status test.
834  outputTest_->reset();
835 
836  // Create the first block in the current Krylov basis (residual).
837  if (R_ == Teuchos::null)
838  R_ = MVT::Clone( *(problem_->getRHS()), 1 );
839 
840  problem_->computeCurrResVec( &*R_ );
841 
842 
843  // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
844  // TODO: ensure hypothesis right here ... I have to think about use cases.
845 
846  if( U_ != Teuchos::null ){
847  // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
848 
849  // possibly over solved equation ... I want residual norms
850  // relative to the initial residual, not what I am about to compute.
851  Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
852  std::vector<MagnitudeType> rnorm0(1);
853  MVT::MvNorm( *R_, rnorm0 ); // rnorm0 = norm(R_);
854 
855  // Z := U_'*R_; xo += U_*Z ;R_ -= C_*Z
856  std::cout << "Solver Manager: dimU_ = " << dimU_ << std::endl;
858 
859  Teuchos::RCP<const MV> Uactive, Cactive;
860  std::vector<int> active_columns( dimU_ );
861  for (int i=0; i < dimU_; ++i) active_columns[i] = i;
862  Uactive = MVT::CloneView(*U_, active_columns);
863  Cactive = MVT::CloneView(*C_, active_columns);
864 
865  if( debug ){
866  std::cout << " Solver Manager : check duality of seed basis " << std::endl;
868  MVT::MvTransMv( one, *Uactive, *Cactive, H );
869  H.print( std::cout );
870  }
871 
872  MVT::MvTransMv( one, *Uactive, *R_, Z );
873  Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
874  MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU ); // UZ
875  MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec ); // xo += tmp;
876  MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU ); // CZ
877  MVT::MvAddMv( -one, *tempU, one, *R_, *R_ ); // R_ -= tmp;
878  std::vector<MagnitudeType> rnorm(1);
879  MVT::MvNorm( *R_, rnorm );
880  if( rnorm[0] < rnorm0[0] * .001 ){ //reorthogonalize
881  MVT::MvTransMv( one, *Uactive, *R_, Z );
882  MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
883  MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec ); // xo += UZ;
884  MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
885  MVT::MvAddMv( -one, *tempU, one, *R_, *R_ ); // R_ -= CZ;
886  }
887  Uactive = Teuchos::null;
888  Cactive = Teuchos::null;
889  tempU = Teuchos::null;
890  }
891  else {
892  dimU_ = 0;
893  }
894 
895 
896  // Set the new state and initialize the solver.
897  PCPGIterState<ScalarType,MV> pcpgState; // fails if R == null.
898 
899  pcpgState.R = R_;
900  if( U_ != Teuchos::null ) pcpgState.U = U_;
901  if( C_ != Teuchos::null ) pcpgState.C = C_;
902  if( dimU_ > 0 ) pcpgState.curDim = dimU_;
903  pcpg_iter->initialize(pcpgState);
904 
905  // treat initialize() exceptions here? how to use try-catch-throw? DMD
906 
907  // Get the current number of deflated blocks with the PCPG iteration
908  dimU_ = pcpgState.curDim;
909  if( !dimU_ ) printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
910  pcpg_iter->resetNumIters();
911 
912  if( dimU_ > savedBlocks_ )
913  std::cout << "Error: dimU_ = " << dimU_ << " > savedBlocks_ = " << savedBlocks_ << std::endl;
914 
915  while(1) { // dummy loop for break
916 
917  // tell pcpg_iter to iterate
918  try {
919  if( debug ) printf("********** Calling iterate...\n");
920  pcpg_iter->iterate();
921 
923  //
924  // check convergence first
925  //
927  if ( convTest_->getStatus() == Passed ) {
928  // we have convergence
929  break; // break from while(1){pcpg_iter->iterate()}
930  }
932  //
933  // check for maximum iterations
934  //
936  else if ( maxIterTest_->getStatus() == Passed ) {
937  // we don't have convergence
938  isConverged = false;
939  break; // break from while(1){pcpg_iter->iterate()}
940  }
941  else {
942 
944  //
945  // we returned from iterate(), but none of our status tests Passed.
946  // Something is wrong, and it is probably the developers fault.
947  //
949 
950  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
951  "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
952  } // end if
953  } // end try
954  catch (const PCPGIterOrthoFailure &e) {
955 
956  // Check to see if the most recent solution yielded convergence.
957  sTest_->checkStatus( &*pcpg_iter );
958  if (convTest_->getStatus() != Passed)
959  isConverged = false;
960  break;
961  }
962  catch (const std::exception &e) {
963  printer_->stream(Errors) << "Error! Caught exception in PCPGIter::iterate() at iteration "
964  << pcpg_iter->getNumIters() << std::endl
965  << e.what() << std::endl;
966  throw;
967  }
968  } // end of while(1)
969 
970  // Update the linear problem.
971  Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
972  problem_->updateSolution( update, true );
973 
974  // Inform the linear problem that we are finished with this block linear system.
975  problem_->setCurrLS();
976 
977  // Get the state. How did pcpgState die?
978  PCPGIterState<ScalarType,MV> oldState = pcpg_iter->getState();
979 
980  dimU_ = oldState.curDim;
981  int q = oldState.prevUdim;
982 
983  std::cout << "SolverManager: dimU_ " << dimU_ << " prevUdim= " << q << std::endl;
984 
985  if( q > deflatedBlocks_ )
986  std::cout << "SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
987 
988  int rank;
989  if( dimU_ > q ){ // Orthogonalize [U;C](:,prevUdim:dimU_)
990  //Given the seed space U and C = A U for some symmetric positive definite A,
991  //find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU
992 
993  //oldState.D->print( std::cout ); D = diag( C'*U )
994 
995  U_ = oldState.U; //MVT::MvPrint( *U_, std::cout );
996  C_ = oldState.C; //MVT::MvPrint( *C_, std::cout );
997  rank = ARRQR(dimU_,q, *oldState.D );
998  if( rank < dimU_ ) {
999  std::cout << " rank decreased in ARRQR, something to do? " << std::endl;
1000  }
1001  dimU_ = rank;
1002 
1003  } // Now U_ and C_ = AU are dual bases.
1004 
1005  if( dimU_ > deflatedBlocks_ ){
1006 
1007  if( !deflatedBlocks_ ){
1008  U_ = Teuchos::null;
1009  C_ = Teuchos::null;
1010  dimU_ = deflatedBlocks_;
1011  break;
1012  }
1013 
1014  bool Harmonic = false; // (Harmonic) Ritz vectors
1015 
1016  Teuchos::RCP<MV> Uorth;
1017 
1018  std::vector<int> active_cols( dimU_ );
1019  for (int i=0; i < dimU_; ++i) active_cols[i] = i;
1020 
1021  if( Harmonic ){
1022  Uorth = MVT::CloneCopy(*C_, active_cols);
1023  }
1024  else{
1025  Uorth = MVT::CloneCopy(*U_, active_cols);
1026  }
1027 
1028  // Explicitly construct Q and R factors
1030  rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,false));
1031  Uorth = Teuchos::null;
1032  // TODO: During the previous solve, the matrix that normalizes U(1:q) was computed and discarded.
1033  // One might save it, reuse it here, and just normalize columns U(q+1:dimU_) here.
1034 
1035  // throw an error if U is both A-orthonormal and rank deficient
1037  "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1038 
1039 
1040  // R VT' = Ur S,
1041  Teuchos::SerialDenseMatrix<int,ScalarType> VT; // Not referenced
1042  Teuchos::SerialDenseMatrix<int,ScalarType> Ur; // Not referenced
1043  int lwork = 5*dimU_; // minimal, extra computation < 67*dimU_
1044  int info = 0; // Hermite
1045  int lrwork = 1;
1046  if( problem_->isHermitian() ) lrwork = dimU_;
1047  std::vector<ScalarType> work(lwork); //
1048  std::vector<ScalarType> Svec(dimU_); //
1049  std::vector<ScalarType> rwork(lrwork);
1050  lapack.GESVD('N', 'O',
1051  R.numRows(),R.numCols(),R.values(), R.numRows(),
1052  &Svec[0],
1053  Ur.values(),1,
1054  VT.values(),1, // Output: VT stored in R
1055  &work[0], lwork,
1056  &rwork[0], &info);
1057 
1059  "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
1060 
1061  if( work[0] != 67. * dimU_ )
1062  std::cout << " SVD " << dimU_ << " lwork " << work[0] << std::endl;
1063  for( int i=0; i< dimU_; i++)
1064  std::cout << i << " " << Svec[i] << std::endl;
1065 
1067 
1068  int startRow = 0, startCol = 0;
1069  if( Harmonic )
1070  startCol = dimU_ - deflatedBlocks_;
1071 
1073  wholeV,
1074  wholeV.numRows(),
1075  deflatedBlocks_,
1076  startRow,
1077  startCol);
1078  std::vector<int> active_columns( dimU_ );
1079  std::vector<int> def_cols( deflatedBlocks_ );
1080  for (int i=0; i < dimU_; ++i) active_columns[i] = i;
1081  for (int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
1082 
1083  Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
1084  Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
1085  MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive ); // U:= U*V
1086  Ucopy = Teuchos::null;
1087  Uactive = Teuchos::null;
1088  Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
1089  Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
1090  MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive ); // C:= C*V
1091  Ccopy = Teuchos::null;
1092  Cactive = Teuchos::null;
1093  dimU_ = deflatedBlocks_;
1094  }
1095  printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << dimU_ << std::endl << std::endl;
1096 
1097  // Inform the linear problem that we are finished with this block linear system.
1098  problem_->setCurrLS();
1099 
1100  // Update indices for the linear systems to be solved.
1101  numRHS2Solve -= 1;
1102  if ( numRHS2Solve > 0 ) {
1103  currIdx[0]++;
1104 
1105  // Set the next indices.
1106  problem_->setLSIndex( currIdx );
1107  }
1108  else {
1109  currIdx.resize( numRHS2Solve );
1110  }
1111  }// while ( numRHS2Solve > 0 )
1112  }
1113 
1114  // print final summary
1115  sTest_->print( printer_->stream(FinalSummary) );
1116 
1117  // print timing information
1118 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1119  // Calling summarize() can be expensive, so don't call unless the
1120  // user wants to print out timing details. summarize() will do all
1121  // the work even if it's passed a "black hole" output stream.
1122  if (verbosity_ & TimingDetails)
1123  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1124 #endif
1125 
1126  // Save the convergence test value ("achieved tolerance") for this solve.
1127  {
1128  using Teuchos::rcp_dynamic_cast;
1129  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1130  // testValues is nonnull and not persistent.
1131  const std::vector<MagnitudeType>* pTestValues =
1132  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1133 
1134  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1135  "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1136  "method returned NULL. Please report this bug to the Belos developers.");
1137 
1138  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1139  "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1140  "method returned a vector of length zero. Please report this bug to the "
1141  "Belos developers.");
1142 
1143  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1144  // achieved tolerances for all vectors in the current solve(), or
1145  // just for the vectors from the last deflation?
1146  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1147  }
1148 
1149  // get iteration information for this solve
1150  numIters_ = maxIterTest_->getNumIters();
1151 
1152  if (!isConverged) {
1153  return Unconverged; // return from PCPGSolMgr::solve()
1154  }
1155  return Converged; // return from PCPGSolMgr::solve()
1156 }
1157 
1158 // A-orthogonalize the Seed Space
1159 // Note that Anasazi::GenOrthoManager provides simplified versions of the algorithm,
1160 // that are not rank revealing, and are not designed for PCPG in other ways too.
1161 template<class ScalarType, class MV, class OP>
1163 {
1164  using Teuchos::RCP;
1165  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1166  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1167 
1168  // Allocate memory for scalars.
1172  std::vector<int> curind(1);
1173  std::vector<int> ipiv(p - q); // RRQR Pivot indices
1174  std::vector<ScalarType> Pivots(p); // RRQR Pivots
1175  int i, imax, j, l;
1176  ScalarType rteps = 1.5e-8;
1177 
1178  // Scale such that diag( U'C) = I
1179  for( i = q ; i < p ; i++ ){
1180  ipiv[i-q] = i;
1181  curind[0] = i;
1182  RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1183  RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1184  anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
1185  MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1186  MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1187  Pivots[i] = one;
1188  }
1189 
1190  for( i = q ; i < p ; i++ ){
1191  if( q < i && i < p-1 ){ // Find the largest pivot
1192  imax = i;
1193  l = ipiv[imax-q];
1194  for( j = i+1 ; j < p ; j++ ){
1195  const int k = ipiv[j-q];
1196  if( Pivots[k] > Pivots[l] ){
1197  imax = j;
1198  l = k;
1199  }
1200  } // end for
1201  if( imax > i ){
1202  l = ipiv[imax-q]; // swap ipiv( imax ) and ipiv(i+1)
1203  ipiv[imax-q] = ipiv[i-q];
1204  ipiv[i-q] = l;
1205  }
1206  } // largest pivot found
1207  int k = ipiv[i-q];
1208 
1209  if( Pivots[k] > 1.5625e-2 ){
1210  anorm(0,0) = Pivots[k]; // A-norm of u
1211  }
1212  else{ // anorm(0,0) = sqrt( U(:,k)'*C(:,k) );
1213  curind[0] = k;
1214  RCP<const MV> P = MVT::CloneView(*U_,curind);
1215  RCP<const MV> AP = MVT::CloneView(*C_,curind);
1216  MVT::MvTransMv( one, *P, *AP, anorm );
1217  anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
1218  }
1219  if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
1220  /*
1221  C(:,k) = A*U(:,k); % Change C
1222  fixC = U(:, ipiv(1:i-1) )'*C(:,k);
1223  U(:,k) = U(:,k) - U(:, ipiv(1:i-1) )*fixC;
1224  C(:,k) = C(:,k) - C(:, ipiv(1:i-1) )*fixC;
1225  anorm = sqrt( U(:,k)'*C(:,k) );
1226  */
1227  std::cout << "ARRQR: Bad case not implemented" << std::endl;
1228  }
1229  if( anorm(0,0) < rteps ){ // rank [U;C] = i-1
1230  std::cout << "ARRQR : deficient case not implemented " << std::endl;
1231  //U = U(:, ipiv(1:i-1) );
1232  //C = C(:, ipiv(1:i-1) );
1233  p = q + i;
1234  // update curDim_ in State
1235  break;
1236  }
1237  curind[0] = k;
1238  RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1239  RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1240  MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P ); // U(:,k) = U(:,k)/anorm;
1241  MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP ); // C(:,k) = C(:,k)/anorm;
1242  P = Teuchos::null;
1243  AP = Teuchos::null;
1244  Pivots[k] = one; // delete, for diagonostic purposes
1245  P = MVT::CloneViewNonConst(*U_,curind); // U(:,k)
1246  AP = MVT::CloneViewNonConst(*C_,curind); // C(:,k)
1247  for( j = i+1 ; j < p ; j++ ){
1248  l = ipiv[j-q]; // ahhh
1249  curind[0] = l;
1250  RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind); // segmentation fault, j=i+1=5
1251  MVT::MvTransMv( one, *Q, *AP, alpha); // alpha(0,0) = U(:,l)'*C(:,k);
1252  MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q ); // U(:,l) -= U(:,k) * alpha(0,0);
1253  Q = Teuchos::null;
1254  RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
1255  MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ ); // C(:,l) -= C(:,l) - C(:,k) * alpha(0,0);
1256  AQ = Teuchos::null;
1257  gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
1258  if( gamma(0,0) > 0){
1259  Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
1260  }
1261  else {
1262  Pivots[l] = zero; //rank deficiency revealed
1263  }
1264  }
1265  }
1266  return p;
1267 }
1268 
1269 // The method returns a string describing the solver manager.
1270 template<class ScalarType, class MV, class OP>
1272 {
1273  std::ostringstream oss;
1274  oss << "Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1275  oss << "{";
1276  oss << "Ortho Type='"<<orthoType_;
1277  oss << "}";
1278  return oss.str();
1279 }
1280 
1281 } // end Belos namespace
1282 
1283 #endif /* BELOS_PCPG_SOLMGR_HPP */
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Collection of types and exceptions used within the Belos solvers.
int prevUdim
Number of block columns in matrices C and U before current iteration.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
PCPGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Class which manages the output and verbosity of the Belos solvers.
static const bool scalarTypeIsSupported
ScalarType * values() const
bool is_null(const boost::shared_ptr< T > &p)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Belos concrete class to iterate Preconditioned Conjugate Projected Gradients.
int numIters_
Number of iterations taken by the last solve() invocation.
PCPGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonormal...
Teuchos::RCP< MV > R
The current residual.
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)
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
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)
int maxIters_
Maximum iteration count (read from parameter list).
virtual void print(std::ostream &os) const
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.
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.
int getNumIters() const
Get the iteration count for the most recent call to solve().
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.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
OperatorTraits< ScalarType, MV, OP > OPT
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:205
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. ...
OrdinalType numRows() const
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.
PCPGSolMgrLAPACKFailure(const std::string &what_arg)
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
static const Teuchos::RCP< std::ostream > outputStream_default_
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:154
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
Teuchos::RCP< OutputManager< ScalarType > > printer_
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::RCP< MV > U
The recycled subspace.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
PCPGIterOrthoFailure is thrown when the PCPGIter object is unable to compute independent direction ve...
Teuchos::ScalarTraits< ScalarType > SCT
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
PCPG iterative linear solver.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
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
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Details::SolverManagerRequiresRealLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
PCPGSolMgrRecyclingFailure(const std::string &what_arg)
OrdinalType numCols() const
This class implements the PCPG iteration, where a single-std::vector Krylov subspace is constructed...
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< Teuchos::ParameterList > params_
PCPGSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::ScalarTraits< MagnitudeType > MT