Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosPCPGSolMgr.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_PCPG_SOLMGR_HPP
11 #define BELOS_PCPG_SOLMGR_HPP
12 
16 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 #include "BelosLinearProblem.hpp"
20 #include "BelosSolverManager.hpp"
21 
22 #include "BelosPCPGIter.hpp"
23 
27 #include "BelosStatusTestCombo.hpp"
29 #include "BelosOutputManager.hpp"
30 #include "Teuchos_LAPACK.hpp"
31 #ifdef BELOS_TEUCHOS_TIME_MONITOR
32 # include "Teuchos_TimeMonitor.hpp"
33 #endif
34 #if defined(HAVE_TEUCHOSCORE_CXX11)
35 # include <type_traits>
36 #endif // defined(HAVE_TEUCHOSCORE_CXX11)
37 
45 namespace Belos {
46 
48 
49 
57  PCPGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
58  {}};
59 
65  class PCPGSolMgrOrthoFailure : public BelosError {public:
66  PCPGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
67  {}};
68 
75  class PCPGSolMgrLAPACKFailure : public BelosError {public:
76  PCPGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
77  {}};
78 
80 
81 
102 
103  // Partial specialization for complex ScalarType.
104  // This contains a trivial implementation.
105  // See discussion in the class documentation above.
106  //
107  // FIXME (mfh 09 Sep 2015) This also is a stub for types other than
108  // float or double.
109  template<class ScalarType, class MV, class OP,
110  const bool supportsScalarType =
113  class PCPGSolMgr :
114  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
115  Belos::Details::LapackSupportsScalar<ScalarType>::value &&
116  ! Teuchos::ScalarTraits<ScalarType>::isComplex>
117  {
118  static const bool scalarTypeIsSupported =
121  typedef Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
122  scalarTypeIsSupported> base_type;
123 
124  public:
126  base_type ()
127  {}
130  base_type ()
131  {}
132  virtual ~PCPGSolMgr () {}
133 
137  }
138  };
139 
140  template<class ScalarType, class MV, class OP>
141  class PCPGSolMgr<ScalarType, MV, OP, true> :
142  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, true> {
143  private:
147  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
149 
150  public:
152 
153 
160  PCPGSolMgr();
161 
199 
201  virtual ~PCPGSolMgr() {};
202 
206  }
208 
210 
211 
215  return *problem_;
216  }
217 
220  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
221 
225 
232  return Teuchos::tuple(timerSolve_);
233  }
234 
240  MagnitudeType achievedTol() const {
241  return achievedTol_;
242  }
243 
245  int getNumIters() const {
246  return numIters_;
247  }
248 
251  bool isLOADetected() const { return false; }
252 
254 
256 
257 
259  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
260 
262  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
263 
265 
267 
268 
272  void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
274 
276 
277 
295  ReturnType solve();
296 
298 
301 
303  std::string description() const;
304 
306 
307  private:
308 
309  // In the A-inner product, perform an RRQR decomposition without using A unless absolutely necessary. Given
310  // 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.
311  int ARRQR(int numVecs, int numOrthVecs, const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
312 
313  // Linear problem.
315 
316  // Output manager.
318  Teuchos::RCP<std::ostream> outputStream_;
319 
320  // Status test.
325 
326  // Orthogonalization manager.
328 
329  // Current parameter list.
331 
332  // Default solver values.
333  static constexpr int maxIters_default_ = 1000;
334  static constexpr int deflatedBlocks_default_ = 2;
335  static constexpr int savedBlocks_default_ = 16;
336  static constexpr int verbosity_default_ = Belos::Errors;
337  static constexpr int outputStyle_default_ = Belos::General;
338  static constexpr int outputFreq_default_ = -1;
339  static constexpr const char * label_default_ = "Belos";
340  static constexpr const char * orthoType_default_ = "ICGS";
341 
342  //
343  // Current solver values.
344  //
345 
347  MagnitudeType convtol_;
348 
350  MagnitudeType orthoKappa_;
351 
353  MagnitudeType achievedTol_;
354 
356  int numIters_;
357 
359  int maxIters_;
360 
361  int deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
362  std::string orthoType_;
363 
364  // Recycled subspace, its image and the residual
365  Teuchos::RCP<MV> U_, C_, R_;
366 
367  // Actual dimension of current recycling subspace (<= savedBlocks_ )
368  int dimU_;
369 
370  // Timers.
371  std::string label_;
372  Teuchos::RCP<Teuchos::Time> timerSolve_;
373 
374  // Internal state variables.
375  bool isSet_;
376  };
377 
378 
379 // Empty Constructor
380 template<class ScalarType, class MV, class OP>
382  outputStream_(Teuchos::rcpFromRef(std::cout)),
383  convtol_(DefaultSolverParameters::convTol),
384  orthoKappa_(DefaultSolverParameters::orthoKappa),
385  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
386  numIters_(0),
387  maxIters_(maxIters_default_),
388  deflatedBlocks_(deflatedBlocks_default_),
389  savedBlocks_(savedBlocks_default_),
390  verbosity_(verbosity_default_),
391  outputStyle_(outputStyle_default_),
392  outputFreq_(outputFreq_default_),
393  orthoType_(orthoType_default_),
394  dimU_(0),
395  label_(label_default_),
396  isSet_(false)
397 {}
398 
399 
400 // Basic Constructor
401 template<class ScalarType, class MV, class OP>
405  problem_(problem),
406  outputStream_(Teuchos::rcpFromRef(std::cout)),
407 
408  convtol_(DefaultSolverParameters::convTol),
409  orthoKappa_(DefaultSolverParameters::orthoKappa),
410  achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
411  numIters_(0),
412  maxIters_(maxIters_default_),
413  deflatedBlocks_(deflatedBlocks_default_),
414  savedBlocks_(savedBlocks_default_),
415  verbosity_(verbosity_default_),
416  outputStyle_(outputStyle_default_),
417  outputFreq_(outputFreq_default_),
418  orthoType_(orthoType_default_),
419  dimU_(0),
420  label_(label_default_),
421  isSet_(false)
422 {
424  problem_.is_null (), std::invalid_argument,
425  "Belos::PCPGSolMgr two-argument constructor: "
426  "'problem' is null. You must supply a non-null Belos::LinearProblem "
427  "instance when calling this constructor.");
428 
429  if (! pl.is_null ()) {
430  // Set the parameters using the list that was passed in.
431  setParameters (pl);
432  }
433 }
434 
435 
436 template<class ScalarType, class MV, class OP>
438 {
439  // Create the internal parameter list if ones doesn't already exist.
440  if (params_ == Teuchos::null) {
441  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
442  }
443  else {
444  params->validateParameters(*getValidParameters());
445  }
446 
447  // Check for maximum number of iterations
448  if (params->isParameter("Maximum Iterations")) {
449  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
450 
451  // Update parameter in our list and in status test.
452  params_->set("Maximum Iterations", maxIters_);
453  if (maxIterTest_!=Teuchos::null)
454  maxIterTest_->setMaxIters( maxIters_ );
455  }
456 
457  // Check for the maximum numbers of saved and deflated blocks.
458  if (params->isParameter("Num Saved Blocks")) {
459  savedBlocks_ = params->get("Num Saved Blocks",savedBlocks_default_);
460  TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
461  "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
462 
463  // savedBlocks > number of matrix rows and columns, not known in parameters.
464  //TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ >= maxIters_, std::invalid_argument,
465  //"Belos::PCPGSolMgr: \"Num Saved Blocks\" must be less than \"Maximum Iterations\".");
466 
467  // Update parameter in our list.
468  params_->set("Num Saved Blocks", savedBlocks_);
469  }
470  if (params->isParameter("Num Deflated Blocks")) {
471  deflatedBlocks_ = params->get("Num Deflated Blocks",deflatedBlocks_default_);
472  TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
473  "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
474 
475  TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
476  "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
477 
478  // Update parameter in our list.
479  // The static_cast is for clang link issues with the constexpr before c++17
480  params_->set("Num Deflated Blocks", static_cast<int>(deflatedBlocks_));
481  }
482 
483  // Check to see if the timer label changed.
484  if (params->isParameter("Timer Label")) {
485  std::string tempLabel = params->get("Timer Label", label_default_);
486 
487  // Update parameter in our list and solver timer
488  if (tempLabel != label_) {
489  label_ = tempLabel;
490  params_->set("Timer Label", label_);
491  std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
492 #ifdef BELOS_TEUCHOS_TIME_MONITOR
493  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
494 #endif
495  if (ortho_ != Teuchos::null) {
496  ortho_->setLabel( label_ );
497  }
498  }
499  }
500 
501  // Check for a change in verbosity level
502  if (params->isParameter("Verbosity")) {
503  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
504  verbosity_ = params->get("Verbosity", verbosity_default_);
505  } else {
506  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
507  }
508 
509  // Update parameter in our list.
510  params_->set("Verbosity", verbosity_);
511  if (printer_ != Teuchos::null)
512  printer_->setVerbosity(verbosity_);
513  }
514 
515  // Check for a change in output style
516  if (params->isParameter("Output Style")) {
517  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
518  outputStyle_ = params->get("Output Style", outputStyle_default_);
519  } else {
520  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
521  }
522 
523  // Reconstruct the convergence test if the explicit residual test is not being used.
524  params_->set("Output Style", outputStyle_);
525  outputTest_ = Teuchos::null;
526  }
527 
528  // output stream
529  if (params->isParameter("Output Stream")) {
530  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
531 
532  // Update parameter in our list.
533  params_->set("Output Stream", outputStream_);
534  if (printer_ != Teuchos::null)
535  printer_->setOStream( outputStream_ );
536  }
537 
538  // frequency level
539  if (verbosity_ & Belos::StatusTestDetails) {
540  if (params->isParameter("Output Frequency")) {
541  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
542  }
543 
544  // Update parameter in out list and output status test.
545  params_->set("Output Frequency", outputFreq_);
546  if (outputTest_ != Teuchos::null)
547  outputTest_->setOutputFrequency( outputFreq_ );
548  }
549 
550  // Create output manager if we need to.
551  if (printer_ == Teuchos::null) {
552  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
553  }
554 
555  // Check if the orthogonalization changed.
556  bool changedOrthoType = false;
557  if (params->isParameter("Orthogonalization")) {
558  std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
559  if (tempOrthoType != orthoType_) {
560  orthoType_ = tempOrthoType;
561  changedOrthoType = true;
562  }
563  }
564  params_->set("Orthogonalization", orthoType_);
565 
566  // Check which orthogonalization constant to use.
567  if (params->isParameter("Orthogonalization Constant")) {
568  if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
569  orthoKappa_ = params->get ("Orthogonalization Constant",
570  static_cast<MagnitudeType> (DefaultSolverParameters::orthoKappa));
571  }
572  else {
573  orthoKappa_ = params->get ("Orthogonalization Constant",
575  }
576 
577  // Update parameter in our list.
578  params_->set("Orthogonalization Constant",orthoKappa_);
579  if (orthoType_=="DGKS") {
580  if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
581  Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
582  }
583  }
584  }
585 
586  // Create orthogonalization manager if we need to.
587  if (ortho_ == Teuchos::null || changedOrthoType) {
589  Teuchos::RCP<Teuchos::ParameterList> paramsOrtho; // can be null
590  if (orthoType_=="DGKS" && orthoKappa_ > 0) {
591  paramsOrtho = Teuchos::rcp(new Teuchos::ParameterList());
592  paramsOrtho->set ("depTol", orthoKappa_ );
593  }
594 
595  ortho_ = factory.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
596  }
597 
598  // Convergence
599  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
600  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
601 
602  // Check for convergence tolerance
603  if (params->isParameter("Convergence Tolerance")) {
604  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
605  convtol_ = params->get ("Convergence Tolerance",
606  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
607  }
608  else {
609  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
610  }
611 
612  // Update parameter in our list and residual tests.
613  params_->set("Convergence Tolerance", convtol_);
614  if (convTest_ != Teuchos::null)
615  convTest_->setTolerance( convtol_ );
616  }
617 
618  // Create status tests if we need to.
619 
620  // Basic test checks maximum iterations and native residual.
621  if (maxIterTest_ == Teuchos::null)
622  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
623 
624  if (convTest_ == Teuchos::null)
625  convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
626 
627  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
628 
629  // Create the status test output class.
630  // This class manages and formats the output from the status test.
631  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
632  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
633 
634  // Set the solver string for the output test
635  std::string solverDesc = " PCPG ";
636  outputTest_->setSolverDesc( solverDesc );
637 
638  // Create the timer if we need to.
639  if (timerSolve_ == Teuchos::null) {
640  std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
641 #ifdef BELOS_TEUCHOS_TIME_MONITOR
642  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
643 #endif
644  }
645 
646  // Inform the solver manager that the current parameters were set.
647  isSet_ = true;
648 }
649 
650 
651 template<class ScalarType, class MV, class OP>
654 {
656  if (is_null(validPL)) {
657  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
658  // Set all the valid parameters and their default values.
659  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
660  "The relative residual tolerance that needs to be achieved by the\n"
661  "iterative solver in order for the linear system to be declared converged.");
662  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
663  "The maximum number of iterations allowed for each\n"
664  "set of RHS solved.");
665  pl->set("Num Deflated Blocks", static_cast<int>(deflatedBlocks_default_),
666  "The maximum number of vectors in the seed subspace." );
667  pl->set("Num Saved Blocks", static_cast<int>(savedBlocks_default_),
668  "The maximum number of vectors saved from old Krylov subspaces." );
669  pl->set("Verbosity", static_cast<int>(verbosity_default_),
670  "What type(s) of solver information should be outputted\n"
671  "to the output stream.");
672  pl->set("Output Style", static_cast<int>(outputStyle_default_),
673  "What style is used for the solver information outputted\n"
674  "to the output stream.");
675  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
676  "How often convergence information should be outputted\n"
677  "to the output stream.");
678  pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
679  "A reference-counted pointer to the output stream where all\n"
680  "solver output is sent.");
681  pl->set("Timer Label", static_cast<const char *>(label_default_),
682  "The string to use as a prefix for the timer labels.");
683  pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
684  "The type of orthogonalization to use: DGKS, ICGS, IMGS");
685  pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(DefaultSolverParameters::orthoKappa),
686  "The constant used by DGKS orthogonalization to determine\n"
687  "whether another step of classical Gram-Schmidt is necessary.");
688  validPL = pl;
689  }
690  return validPL;
691 }
692 
693 
694 // solve()
695 template<class ScalarType, class MV, class OP>
697 
698  // Set the current parameters if are not set already.
699  if (!isSet_) { setParameters( params_ ); }
700 
702  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
703  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
704 
706  "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
707 
709  "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
710 
711  // Create indices for the linear systems to be solved.
712  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
713  std::vector<int> currIdx(1);
714  currIdx[0] = 0;
715 
716  bool debug = false;
717 
718  // Inform the linear problem of the current linear system to solve.
719  problem_->setLSIndex( currIdx ); // block size == 1
720 
721  // Assume convergence is achieved, then let any failed convergence set this to false.
722  bool isConverged = true;
723 
725  // PCPG iteration parameter list
727  plist.set("Saved Blocks", savedBlocks_);
728  plist.set("Block Size", 1);
729  plist.set("Keep Diagonal", true);
730  plist.set("Initialize Diagonal", true);
731 
733  // PCPG solver
734 
736  pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
737  // Number of iterations required to generate initial recycle space (if needed)
738 
739  // Enter solve() iterations
740  {
741 #ifdef BELOS_TEUCHOS_TIME_MONITOR
742  Teuchos::TimeMonitor slvtimer(*timerSolve_);
743 #endif
744  while ( numRHS2Solve > 0 ) { // test for quick return
745 
746  // Reset the status test.
747  outputTest_->reset();
748 
749  // Create the first block in the current Krylov basis (residual).
750  if (R_ == Teuchos::null)
751  R_ = MVT::Clone( *(problem_->getRHS()), 1 );
752 
753  problem_->computeCurrResVec( &*R_ );
754 
755 
756  // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
757  // TODO: ensure hypothesis right here ... I have to think about use cases.
758 
759  if( U_ != Teuchos::null ){
760  // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
761 
762  // possibly over solved equation ... I want residual norms
763  // relative to the initial residual, not what I am about to compute.
764  Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
765  std::vector<MagnitudeType> rnorm0(1);
766  MVT::MvNorm( *R_, rnorm0 ); // rnorm0 = norm(R_);
767 
768  // Z := U_'*R_; xo += U_*Z ;R_ -= C_*Z
769  std::cout << "Solver Manager: dimU_ = " << dimU_ << std::endl;
771 
772  Teuchos::RCP<const MV> Uactive, Cactive;
773  std::vector<int> active_columns( dimU_ );
774  for (int i=0; i < dimU_; ++i) active_columns[i] = i;
775  Uactive = MVT::CloneView(*U_, active_columns);
776  Cactive = MVT::CloneView(*C_, active_columns);
777 
778  if( debug ){
779  std::cout << " Solver Manager : check duality of seed basis " << std::endl;
781  MVT::MvTransMv( one, *Uactive, *Cactive, H );
782  H.print( std::cout );
783  }
784 
785  MVT::MvTransMv( one, *Uactive, *R_, Z );
786  Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
787  MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU ); // UZ
788  MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec ); // xo += tmp;
789  MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU ); // CZ
790  MVT::MvAddMv( -one, *tempU, one, *R_, *R_ ); // R_ -= tmp;
791  std::vector<MagnitudeType> rnorm(1);
792  MVT::MvNorm( *R_, rnorm );
793  if( rnorm[0] < rnorm0[0] * .001 ){ //reorthogonalize
794  MVT::MvTransMv( one, *Uactive, *R_, Z );
795  MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
796  MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec ); // xo += UZ;
797  MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
798  MVT::MvAddMv( -one, *tempU, one, *R_, *R_ ); // R_ -= CZ;
799  }
800  Uactive = Teuchos::null;
801  Cactive = Teuchos::null;
802  tempU = Teuchos::null;
803  }
804  else {
805  dimU_ = 0;
806  }
807 
808 
809  // Set the new state and initialize the solver.
810  PCPGIterState<ScalarType,MV> pcpgState; // fails if R == null.
811 
812  pcpgState.R = R_;
813  if( U_ != Teuchos::null ) pcpgState.U = U_;
814  if( C_ != Teuchos::null ) pcpgState.C = C_;
815  if( dimU_ > 0 ) pcpgState.curDim = dimU_;
816  pcpg_iter->initialize(pcpgState);
817 
818  // treat initialize() exceptions here? how to use try-catch-throw? DMD
819 
820  // Get the current number of deflated blocks with the PCPG iteration
821  dimU_ = pcpgState.curDim;
822  if( !dimU_ ) printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
823  pcpg_iter->resetNumIters();
824 
825  if( dimU_ > savedBlocks_ )
826  std::cout << "Error: dimU_ = " << dimU_ << " > savedBlocks_ = " << savedBlocks_ << std::endl;
827 
828  while(1) { // dummy loop for break
829 
830  // tell pcpg_iter to iterate
831  try {
832  if( debug ) printf("********** Calling iterate...\n");
833  pcpg_iter->iterate();
834 
836  //
837  // check convergence first
838  //
840  if ( convTest_->getStatus() == Passed ) {
841  // we have convergence
842  break; // break from while(1){pcpg_iter->iterate()}
843  }
845  //
846  // check for maximum iterations
847  //
849  else if ( maxIterTest_->getStatus() == Passed ) {
850  // we don't have convergence
851  isConverged = false;
852  break; // break from while(1){pcpg_iter->iterate()}
853  }
854  else {
855 
857  //
858  // we returned from iterate(), but none of our status tests Passed.
859  // Something is wrong, and it is probably the developers fault.
860  //
862 
863  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
864  "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
865  } // end if
866  } // end try
867  catch (const StatusTestNaNError& e) {
868  // A NaN was detected in the solver. Set the solution to zero and return unconverged.
869  achievedTol_ = MT::one();
870  Teuchos::RCP<MV> X = problem_->getLHS();
871  MVT::MvInit( *X, SCT::zero() );
872  printer_->stream(Warnings) << "Belos::PCPG::solve(): Warning! NaN has been detected!"
873  << std::endl;
874  return Unconverged;
875  }
876  catch (const std::exception &e) {
877  printer_->stream(Errors) << "Error! Caught exception in PCPGIter::iterate() at iteration "
878  << pcpg_iter->getNumIters() << std::endl
879  << e.what() << std::endl;
880  throw;
881  }
882  } // end of while(1)
883 
884  // Update the linear problem.
885  Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
886  problem_->updateSolution( update, true );
887 
888  // Inform the linear problem that we are finished with this block linear system.
889  problem_->setCurrLS();
890 
891  // Get the state. How did pcpgState die?
892  PCPGIterState<ScalarType,MV> oldState = pcpg_iter->getState();
893 
894  dimU_ = oldState.curDim;
895  int q = oldState.prevUdim;
896 
897  std::cout << "SolverManager: dimU_ " << dimU_ << " prevUdim= " << q << std::endl;
898 
899  if( q > deflatedBlocks_ )
900  std::cout << "SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
901 
902  int rank;
903  if( dimU_ > q ){ // Orthogonalize [U;C](:,prevUdim:dimU_)
904  //Given the seed space U and C = A U for some symmetric positive definite A,
905  //find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU
906 
907  //oldState.D->print( std::cout ); D = diag( C'*U )
908 
909  U_ = oldState.U; //MVT::MvPrint( *U_, std::cout );
910  C_ = oldState.C; //MVT::MvPrint( *C_, std::cout );
911  rank = ARRQR(dimU_,q, *oldState.D );
912  if( rank < dimU_ ) {
913  std::cout << " rank decreased in ARRQR, something to do? " << std::endl;
914  }
915  dimU_ = rank;
916 
917  } // Now U_ and C_ = AU are dual bases.
918 
919  if( dimU_ > deflatedBlocks_ ){
920 
921  if( !deflatedBlocks_ ){
922  U_ = Teuchos::null;
923  C_ = Teuchos::null;
924  dimU_ = deflatedBlocks_;
925  break;
926  }
927 
928  bool Harmonic = false; // (Harmonic) Ritz vectors
929 
930  Teuchos::RCP<MV> Uorth;
931 
932  std::vector<int> active_cols( dimU_ );
933  for (int i=0; i < dimU_; ++i) active_cols[i] = i;
934 
935  if( Harmonic ){
936  Uorth = MVT::CloneCopy(*C_, active_cols);
937  }
938  else{
939  Uorth = MVT::CloneCopy(*U_, active_cols);
940  }
941 
942  // Explicitly construct Q and R factors
944  rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,false));
945  Uorth = Teuchos::null;
946  // TODO: During the previous solve, the matrix that normalizes U(1:q) was computed and discarded.
947  // One might save it, reuse it here, and just normalize columns U(q+1:dimU_) here.
948 
949  // throw an error if U is both A-orthonormal and rank deficient
951  "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
952 
953 
954  // R VT' = Ur S,
955  Teuchos::SerialDenseMatrix<int,ScalarType> VT; // Not referenced
956  Teuchos::SerialDenseMatrix<int,ScalarType> Ur; // Not referenced
957  int lwork = 5*dimU_; // minimal, extra computation < 67*dimU_
958  int info = 0; // Hermite
959  int lrwork = 1;
960  if( problem_->isHermitian() ) lrwork = dimU_;
961  std::vector<ScalarType> work(lwork); //
962  std::vector<ScalarType> Svec(dimU_); //
963  std::vector<ScalarType> rwork(lrwork);
964  lapack.GESVD('N', 'O',
965  R.numRows(),R.numCols(),R.values(), R.numRows(),
966  &Svec[0],
967  Ur.values(),1,
968  VT.values(),1, // Output: VT stored in R
969  &work[0], lwork,
970  &rwork[0], &info);
971 
973  "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
974 
975  if( work[0] != 67. * dimU_ )
976  std::cout << " SVD " << dimU_ << " lwork " << work[0] << std::endl;
977  for( int i=0; i< dimU_; i++)
978  std::cout << i << " " << Svec[i] << std::endl;
979 
981 
982  int startRow = 0, startCol = 0;
983  if( Harmonic )
984  startCol = dimU_ - deflatedBlocks_;
985 
987  wholeV,
988  wholeV.numRows(),
989  deflatedBlocks_,
990  startRow,
991  startCol);
992  std::vector<int> active_columns( dimU_ );
993  std::vector<int> def_cols( deflatedBlocks_ );
994  for (int i=0; i < dimU_; ++i) active_columns[i] = i;
995  for (int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
996 
997  Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
998  Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
999  MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive ); // U:= U*V
1000  Ucopy = Teuchos::null;
1001  Uactive = Teuchos::null;
1002  Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
1003  Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
1004  MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive ); // C:= C*V
1005  Ccopy = Teuchos::null;
1006  Cactive = Teuchos::null;
1007  dimU_ = deflatedBlocks_;
1008  }
1009  printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << dimU_ << std::endl << std::endl;
1010 
1011  // Inform the linear problem that we are finished with this block linear system.
1012  problem_->setCurrLS();
1013 
1014  // Update indices for the linear systems to be solved.
1015  numRHS2Solve -= 1;
1016  if ( numRHS2Solve > 0 ) {
1017  currIdx[0]++;
1018 
1019  // Set the next indices.
1020  problem_->setLSIndex( currIdx );
1021  }
1022  else {
1023  currIdx.resize( numRHS2Solve );
1024  }
1025  }// while ( numRHS2Solve > 0 )
1026  }
1027 
1028  // print final summary
1029  sTest_->print( printer_->stream(FinalSummary) );
1030 
1031  // print timing information
1032 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1033  // Calling summarize() can be expensive, so don't call unless the
1034  // user wants to print out timing details. summarize() will do all
1035  // the work even if it's passed a "black hole" output stream.
1036  if (verbosity_ & TimingDetails)
1037  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1038 #endif
1039 
1040  // Save the convergence test value ("achieved tolerance") for this solve.
1041  {
1042  using Teuchos::rcp_dynamic_cast;
1043  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1044  // testValues is nonnull and not persistent.
1045  const std::vector<MagnitudeType>* pTestValues =
1046  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1047 
1048  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1049  "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1050  "method returned NULL. Please report this bug to the Belos developers.");
1051 
1052  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1053  "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
1054  "method returned a vector of length zero. Please report this bug to the "
1055  "Belos developers.");
1056 
1057  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1058  // achieved tolerances for all vectors in the current solve(), or
1059  // just for the vectors from the last deflation?
1060  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1061  }
1062 
1063  // get iteration information for this solve
1064  numIters_ = maxIterTest_->getNumIters();
1065 
1066  if (!isConverged) {
1067  return Unconverged; // return from PCPGSolMgr::solve()
1068  }
1069  return Converged; // return from PCPGSolMgr::solve()
1070 }
1071 
1072 // A-orthogonalize the Seed Space
1073 // Note that Anasazi::GenOrthoManager provides simplified versions of the algorithm,
1074 // that are not rank revealing, and are not designed for PCPG in other ways too.
1075 template<class ScalarType, class MV, class OP>
1077 {
1078  using Teuchos::RCP;
1079  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1080  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1081 
1082  // Allocate memory for scalars.
1086  std::vector<int> curind(1);
1087  std::vector<int> ipiv(p - q); // RRQR Pivot indices
1088  std::vector<ScalarType> Pivots(p); // RRQR Pivots
1089  int i, imax, j, l;
1090  ScalarType rteps = 1.5e-8;
1091 
1092  // Scale such that diag( U'C) = I
1093  for( i = q ; i < p ; i++ ){
1094  ipiv[i-q] = i;
1095  curind[0] = i;
1096  RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1097  RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1098  anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
1099  MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
1100  MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
1101  Pivots[i] = one;
1102  }
1103 
1104  for( i = q ; i < p ; i++ ){
1105  if( q < i && i < p-1 ){ // Find the largest pivot
1106  imax = i;
1107  l = ipiv[imax-q];
1108  for( j = i+1 ; j < p ; j++ ){
1109  const int k = ipiv[j-q];
1110  if( Pivots[k] > Pivots[l] ){
1111  imax = j;
1112  l = k;
1113  }
1114  } // end for
1115  if( imax > i ){
1116  l = ipiv[imax-q]; // swap ipiv( imax ) and ipiv(i+1)
1117  ipiv[imax-q] = ipiv[i-q];
1118  ipiv[i-q] = l;
1119  }
1120  } // largest pivot found
1121  int k = ipiv[i-q];
1122 
1123  if( Pivots[k] > 1.5625e-2 ){
1124  anorm(0,0) = Pivots[k]; // A-norm of u
1125  }
1126  else{ // anorm(0,0) = sqrt( U(:,k)'*C(:,k) );
1127  curind[0] = k;
1128  RCP<const MV> P = MVT::CloneView(*U_,curind);
1129  RCP<const MV> AP = MVT::CloneView(*C_,curind);
1130  MVT::MvTransMv( one, *P, *AP, anorm );
1131  anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
1132  }
1133  if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
1134  /*
1135  C(:,k) = A*U(:,k); % Change C
1136  fixC = U(:, ipiv(1:i-1) )'*C(:,k);
1137  U(:,k) = U(:,k) - U(:, ipiv(1:i-1) )*fixC;
1138  C(:,k) = C(:,k) - C(:, ipiv(1:i-1) )*fixC;
1139  anorm = sqrt( U(:,k)'*C(:,k) );
1140  */
1141  std::cout << "ARRQR: Bad case not implemented" << std::endl;
1142  }
1143  if( anorm(0,0) < rteps ){ // rank [U;C] = i-1
1144  std::cout << "ARRQR : deficient case not implemented " << std::endl;
1145  //U = U(:, ipiv(1:i-1) );
1146  //C = C(:, ipiv(1:i-1) );
1147  p = q + i;
1148  // update curDim_ in State
1149  break;
1150  }
1151  curind[0] = k;
1152  RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
1153  RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
1154  MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P ); // U(:,k) = U(:,k)/anorm;
1155  MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP ); // C(:,k) = C(:,k)/anorm;
1156  P = Teuchos::null;
1157  AP = Teuchos::null;
1158  Pivots[k] = one; // delete, for diagonostic purposes
1159  P = MVT::CloneViewNonConst(*U_,curind); // U(:,k)
1160  AP = MVT::CloneViewNonConst(*C_,curind); // C(:,k)
1161  for( j = i+1 ; j < p ; j++ ){
1162  l = ipiv[j-q]; // ahhh
1163  curind[0] = l;
1164  RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind); // segmentation fault, j=i+1=5
1165  MVT::MvTransMv( one, *Q, *AP, alpha); // alpha(0,0) = U(:,l)'*C(:,k);
1166  MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q ); // U(:,l) -= U(:,k) * alpha(0,0);
1167  Q = Teuchos::null;
1168  RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
1169  MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ ); // C(:,l) -= C(:,l) - C(:,k) * alpha(0,0);
1170  AQ = Teuchos::null;
1171  gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
1172  if( gamma(0,0) > 0){
1173  Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
1174  }
1175  else {
1176  Pivots[l] = zero; //rank deficiency revealed
1177  }
1178  }
1179  }
1180  return p;
1181 }
1182 
1183 // The method returns a string describing the solver manager.
1184 template<class ScalarType, class MV, class OP>
1186 {
1187  std::ostringstream oss;
1188  oss << "Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1189  oss << "{";
1190  oss << "Ortho Type='"<<orthoType_;
1191  oss << "}";
1192  return oss.str();
1193 }
1194 
1195 } // end Belos namespace
1196 
1197 #endif /* BELOS_PCPG_SOLMGR_HPP */
static const double orthoKappa
DGKS orthogonalization constant.
Definition: BelosTypes.hpp:267
ScalarType * values() const
Collection of types and exceptions used within the Belos solvers.
int prevUdim
Number of block columns in matrices C and U before current iteration.
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
PCPGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
Belos concrete class to iterate Preconditioned Conjugate Projected Gradients.
PCPGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonormal...
Teuchos::RCP< MV > R
The current residual.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
static T squareroot(T x)
Teuchos::RCP< MV > C
C = AU, U spans recycled subspace.
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
#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)
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
Structure to contain pointers to PCPGIter state variables.
PCPGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:261
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.
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:174
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
A linear system to solve, and its associated information.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
Class which describes the linear problem to be solved by the iterative solver.
PCPGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
virtual Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const
clone for Inverted Injection (DII)
PCPGSolMgrLAPACKFailure(const std::string &what_arg)
int getNumIters() const
Get the iteration count for the most recent call to solve().
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:123
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
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::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
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
PCPG iterative linear solver.
bool isType(const std::string &name) const
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
Default parameters common to most Belos solvers.
Definition: BelosTypes.hpp:251
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...
PCPGSolMgrLinearProblemFailure(const std::string &what_arg)
OrdinalType numRows() const
bool is_null() const
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.

Generated on Tue Jul 23 2024 09:24:22 for Belos by doxygen 1.8.5