Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosRCGSolMgr.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_RCG_SOLMGR_HPP
11 #define BELOS_RCG_SOLMGR_HPP
12 
17 #include "BelosConfigDefs.hpp"
18 #include "BelosTypes.hpp"
19 
20 #include "BelosLinearProblem.hpp"
21 #include "BelosSolverManager.hpp"
22 
23 #include "BelosRCGIter.hpp"
26 #include "BelosStatusTestCombo.hpp"
28 #include "BelosOutputManager.hpp"
29 #include "Teuchos_LAPACK.hpp"
30 #include "Teuchos_as.hpp"
31 #ifdef BELOS_TEUCHOS_TIME_MONITOR
32 #include "Teuchos_TimeMonitor.hpp"
33 #endif
34 
81 namespace Belos {
82 
84 
85 
92  class RCGSolMgrLinearProblemFailure : public BelosError {public:
93  RCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
94  {}};
95 
102  class RCGSolMgrLAPACKFailure : public BelosError {public:
103  RCGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
104  {}};
105 
107 
108 
109  // Partial specialization for unsupported ScalarType types.
110  // This contains a stub implementation.
111  template<class ScalarType, class MV, class OP,
112  const bool supportsScalarType =
115  class RCGSolMgr :
116  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
117  Belos::Details::LapackSupportsScalar<ScalarType>::value &&
118  ! Teuchos::ScalarTraits<ScalarType>::isComplex>
119  {
120  static const bool scalarTypeIsSupported =
123  typedef Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
124  scalarTypeIsSupported> base_type;
125 
126  public:
128  base_type ()
129  {}
132  base_type ()
133  {}
134  virtual ~RCGSolMgr () {}
135 
139  }
140  };
141 
142  // Partial specialization for real ScalarType.
143  // This contains the actual working implementation of RCG.
144  // See discussion in the class documentation above.
145  template<class ScalarType, class MV, class OP>
146  class RCGSolMgr<ScalarType, MV, OP, true> :
147  public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, true> {
148  private:
152  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
154 
155  public:
156 
158 
159 
165  RCGSolMgr();
166 
190 
192  virtual ~RCGSolMgr() {};
193 
197  }
199 
201 
202 
203  const LinearProblem<ScalarType,MV,OP>& getProblem() const override {
204  return *problem_;
205  }
206 
208  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
209 
212 
219  return Teuchos::tuple(timerSolve_);
220  }
221 
226  MagnitudeType achievedTol() const override {
227  return achievedTol_;
228  }
229 
231  int getNumIters() const override {
232  return numIters_;
233  }
234 
236  bool isLOADetected() const override { return false; }
237 
239 
241 
242 
244  void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override { problem_ = problem; }
245 
247  void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
248 
250 
252 
253 
259  void reset( const ResetType type ) override {
260  if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
261  else if (type & Belos::RecycleSubspace) existU_ = false;
262  }
264 
266 
267 
285  ReturnType solve() override;
286 
288 
291 
293  std::string description() const override;
294 
296 
297  private:
298 
299  // Called by all constructors; Contains init instructions common to all constructors
300  void init();
301 
302  // Computes harmonic eigenpairs of projected matrix created during one cycle.
303  // Y contains the harmonic Ritz vectors corresponding to the recycleBlocks eigenvalues of smallest magnitude.
304  void getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType> &F,
307 
308  // Sort list of n floating-point numbers and return permutation vector
309  void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
310 
311  // Initialize solver state storage
312  void initializeStateStorage();
313 
314  // Linear problem.
316 
317  // Output manager.
319  Teuchos::RCP<std::ostream> outputStream_;
320 
321  // Status test.
326 
327  // Current parameter list.
329 
330  // Default solver values.
331  static constexpr int maxIters_default_ = 1000;
332  static constexpr int blockSize_default_ = 1;
333  static constexpr int numBlocks_default_ = 25;
334  static constexpr int recycleBlocks_default_ = 3;
335  static constexpr bool showMaxResNormOnly_default_ = false;
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 
341  //
342  // Current solver values.
343  //
344 
346  MagnitudeType convtol_;
347 
352  MagnitudeType achievedTol_;
353 
355  int maxIters_;
356 
358  int numIters_;
359 
360  int numBlocks_, recycleBlocks_;
361  bool showMaxResNormOnly_;
362  int verbosity_, outputStyle_, outputFreq_;
363 
365  // Solver State Storage
367  // Search vectors
368  Teuchos::RCP<MV> P_;
369  //
370  // A times current search direction
371  Teuchos::RCP<MV> Ap_;
372  //
373  // Residual vector
374  Teuchos::RCP<MV> r_;
375  //
376  // Preconditioned residual
377  Teuchos::RCP<MV> z_;
378  //
379  // Flag indicating that the recycle space should be used
380  bool existU_;
381  //
382  // Flag indicating that the updated recycle space has been created
383  bool existU1_;
384  //
385  // Recycled subspace and its image
386  Teuchos::RCP<MV> U_, AU_;
387  //
388  // Recycled subspace for next system and its image
389  Teuchos::RCP<MV> U1_;
390  //
391  // Coefficients arising in RCG iteration
395  //
396  // Solutions to local least-squares problems
398  //
399  // The matrix U^T A U
401  //
402  // LU factorization of U^T A U
404  //
405  // Data from LU factorization of UTAU
407  //
408  // The matrix (AU)^T AU
410  //
411  // The scalar r'*z
413  //
414  // Matrices needed for calculation of harmonic Ritz eigenproblem
416  //
417  // Matrices needed for updating recycle space
418  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_;
419  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_;
422  Teuchos::RCP<MV> U1Y1_, PY2_;
424  ScalarType dold;
426 
427  // Timers.
428  std::string label_;
429  Teuchos::RCP<Teuchos::Time> timerSolve_;
430 
431  // Internal state variables.
432  bool params_Set_;
433  };
434 
435 
436 // Empty Constructor
437 template<class ScalarType, class MV, class OP>
439  achievedTol_(0.0),
440  numIters_(0)
441 {
442  init();
443 }
444 
445 // Basic Constructor
446 template<class ScalarType, class MV, class OP>
450  problem_(problem),
451  achievedTol_(0.0),
452  numIters_(0)
453 {
454  init();
455  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
456 
457  // If the parameter list pointer is null, then set the current parameters to the default parameter list.
458  if ( !is_null(pl) ) {
459  setParameters( pl );
460  }
461 }
462 
463 // Common instructions executed in all constructors
464 template<class ScalarType, class MV, class OP>
466 {
467  outputStream_ = Teuchos::rcpFromRef(std::cout);
469  maxIters_ = maxIters_default_;
470  numBlocks_ = numBlocks_default_;
471  recycleBlocks_ = recycleBlocks_default_;
472  verbosity_ = verbosity_default_;
473  outputStyle_= outputStyle_default_;
474  outputFreq_= outputFreq_default_;
475  showMaxResNormOnly_ = showMaxResNormOnly_default_;
476  label_ = label_default_;
477  params_Set_ = false;
478  P_ = Teuchos::null;
479  Ap_ = Teuchos::null;
480  r_ = Teuchos::null;
481  z_ = Teuchos::null;
482  existU_ = false;
483  existU1_ = false;
484  U_ = Teuchos::null;
485  AU_ = Teuchos::null;
486  U1_ = Teuchos::null;
487  Alpha_ = Teuchos::null;
488  Beta_ = Teuchos::null;
489  D_ = Teuchos::null;
490  Delta_ = Teuchos::null;
491  UTAU_ = Teuchos::null;
492  LUUTAU_ = Teuchos::null;
493  ipiv_ = Teuchos::null;
494  AUTAU_ = Teuchos::null;
495  rTz_old_ = Teuchos::null;
496  F_ = Teuchos::null;
497  G_ = Teuchos::null;
498  Y_ = Teuchos::null;
499  L2_ = Teuchos::null;
500  DeltaL2_ = Teuchos::null;
501  AU1TUDeltaL2_ = Teuchos::null;
502  AU1TAU1_ = Teuchos::null;
503  AU1TU1_ = Teuchos::null;
504  AU1TAP_ = Teuchos::null;
505  FY_ = Teuchos::null;
506  GY_ = Teuchos::null;
507  APTAP_ = Teuchos::null;
508  U1Y1_ = Teuchos::null;
509  PY2_ = Teuchos::null;
510  AUTAP_ = Teuchos::null;
511  AU1TU_ = Teuchos::null;
512  dold = 0.;
513 }
514 
515 template<class ScalarType, class MV, class OP>
517 {
518  // Create the internal parameter list if ones doesn't already exist.
519  if (params_ == Teuchos::null) {
520  params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
521  }
522  else {
523  params->validateParameters(*getValidParameters());
524  }
525 
526  // Check for maximum number of iterations
527  if (params->isParameter("Maximum Iterations")) {
528  maxIters_ = params->get("Maximum Iterations",maxIters_default_);
529 
530  // Update parameter in our list and in status test.
531  params_->set("Maximum Iterations", maxIters_);
532  if (maxIterTest_!=Teuchos::null)
533  maxIterTest_->setMaxIters( maxIters_ );
534  }
535 
536  // Check for the maximum number of blocks.
537  if (params->isParameter("Num Blocks")) {
538  numBlocks_ = params->get("Num Blocks",numBlocks_default_);
539  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
540  "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
541 
542  // Update parameter in our list.
543  params_->set("Num Blocks", numBlocks_);
544  }
545 
546  // Check for the maximum number of blocks.
547  if (params->isParameter("Num Recycled Blocks")) {
548  recycleBlocks_ = params->get("Num Recycled Blocks",recycleBlocks_default_);
549  TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
550  "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
551 
552  TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
553  "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
554 
555  // Update parameter in our list.
556  params_->set("Num Recycled Blocks", recycleBlocks_);
557  }
558 
559  // Check to see if the timer label changed.
560  if (params->isParameter("Timer Label")) {
561  std::string tempLabel = params->get("Timer Label", label_default_);
562 
563  // Update parameter in our list and solver timer
564  if (tempLabel != label_) {
565  label_ = tempLabel;
566  params_->set("Timer Label", label_);
567  std::string solveLabel = label_ + ": RCGSolMgr total solve time";
568 #ifdef BELOS_TEUCHOS_TIME_MONITOR
569  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
570 #endif
571  }
572  }
573 
574  // Check for a change in verbosity level
575  if (params->isParameter("Verbosity")) {
576  if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
577  verbosity_ = params->get("Verbosity", verbosity_default_);
578  } else {
579  verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
580  }
581 
582  // Update parameter in our list.
583  params_->set("Verbosity", verbosity_);
584  if (printer_ != Teuchos::null)
585  printer_->setVerbosity(verbosity_);
586  }
587 
588  // Check for a change in output style
589  if (params->isParameter("Output Style")) {
590  if (Teuchos::isParameterType<int>(*params,"Output Style")) {
591  outputStyle_ = params->get("Output Style", outputStyle_default_);
592  } else {
593  outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
594  }
595 
596  // Reconstruct the convergence test if the explicit residual test is not being used.
597  params_->set("Output Style", outputStyle_);
598  outputTest_ = Teuchos::null;
599  }
600 
601  // output stream
602  if (params->isParameter("Output Stream")) {
603  outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
604 
605  // Update parameter in our list.
606  params_->set("Output Stream", outputStream_);
607  if (printer_ != Teuchos::null)
608  printer_->setOStream( outputStream_ );
609  }
610 
611  // frequency level
612  if (verbosity_ & Belos::StatusTestDetails) {
613  if (params->isParameter("Output Frequency")) {
614  outputFreq_ = params->get("Output Frequency", outputFreq_default_);
615  }
616 
617  // Update parameter in out list and output status test.
618  params_->set("Output Frequency", outputFreq_);
619  if (outputTest_ != Teuchos::null)
620  outputTest_->setOutputFrequency( outputFreq_ );
621  }
622 
623  // Create output manager if we need to.
624  if (printer_ == Teuchos::null) {
625  printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
626  }
627 
628  // Convergence
629  typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
630  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
631 
632  // Check for convergence tolerance
633  if (params->isParameter("Convergence Tolerance")) {
634  if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
635  convtol_ = params->get ("Convergence Tolerance",
636  static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
637  }
638  else {
639  convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
640  }
641 
642  // Update parameter in our list and residual tests.
643  params_->set("Convergence Tolerance", convtol_);
644  if (convTest_ != Teuchos::null)
645  convTest_->setTolerance( convtol_ );
646  }
647 
648  if (params->isParameter("Show Maximum Residual Norm Only")) {
649  showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
650 
651  // Update parameter in our list and residual tests
652  params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
653  if (convTest_ != Teuchos::null)
654  convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
655  }
656 
657  // Create status tests if we need to.
658 
659  // Basic test checks maximum iterations and native residual.
660  if (maxIterTest_ == Teuchos::null)
661  maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
662 
663  // Implicit residual test, using the native residual to determine if convergence was achieved.
664  if (convTest_ == Teuchos::null)
665  convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
666 
667  if (sTest_ == Teuchos::null)
668  sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
669 
670  if (outputTest_ == Teuchos::null) {
671 
672  // Create the status test output class.
673  // This class manages and formats the output from the status test.
674  StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
675  outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
676 
677  // Set the solver string for the output test
678  std::string solverDesc = " Recycling CG ";
679  outputTest_->setSolverDesc( solverDesc );
680  }
681 
682  // Create the timer if we need to.
683  if (timerSolve_ == Teuchos::null) {
684  std::string solveLabel = label_ + ": RCGSolMgr total solve time";
685 #ifdef BELOS_TEUCHOS_TIME_MONITOR
686  timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
687 #endif
688  }
689 
690  // Inform the solver manager that the current parameters were set.
691  params_Set_ = true;
692 }
693 
694 
695 template<class ScalarType, class MV, class OP>
698 {
700 
701  // Set all the valid parameters and their default values.
702  if(is_null(validPL)) {
703  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
704  pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
705  "The relative residual tolerance that needs to be achieved by the\n"
706  "iterative solver in order for the linear system to be declared converged.");
707  pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
708  "The maximum number of block iterations allowed for each\n"
709  "set of RHS solved.");
710  pl->set("Block Size", static_cast<int>(blockSize_default_),
711  "Block Size Parameter -- currently must be 1 for RCG");
712  pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
713  "The length of a cycle (and this max number of search vectors kept)\n");
714  pl->set("Num Recycled Blocks", static_cast<int>(recycleBlocks_default_),
715  "The number of vectors in the recycle subspace.");
716  pl->set("Verbosity", static_cast<int>(verbosity_default_),
717  "What type(s) of solver information should be outputted\n"
718  "to the output stream.");
719  pl->set("Output Style", static_cast<int>(outputStyle_default_),
720  "What style is used for the solver information outputted\n"
721  "to the output stream.");
722  pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
723  "How often convergence information should be outputted\n"
724  "to the output stream.");
725  pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
726  "A reference-counted pointer to the output stream where all\n"
727  "solver output is sent.");
728  pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
729  "When convergence information is printed, only show the maximum\n"
730  "relative residual norm when the block size is greater than one.");
731  pl->set("Timer Label", static_cast<const char *>(label_default_),
732  "The string to use as a prefix for the timer labels.");
733  validPL = pl;
734  }
735  return validPL;
736 }
737 
738 // initializeStateStorage
739 template<class ScalarType, class MV, class OP>
741 
742  // Check if there is any multivector to clone from.
743  Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
744  if (rhsMV == Teuchos::null) {
745  // Nothing to do
746  return;
747  }
748  else {
749 
750  // Initialize the state storage
751  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
752  "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
753 
754  // If the subspace has not been initialized before, generate it using the RHS from lp_.
755  if (P_ == Teuchos::null) {
756  P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
757  }
758  else {
759  // Generate P_ by cloning itself ONLY if more space is needed.
760  if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
761  Teuchos::RCP<const MV> tmp = P_;
762  P_ = MVT::Clone( *tmp, numBlocks_+2 );
763  }
764  }
765 
766  // Generate Ap_ only if it doesn't exist
767  if (Ap_ == Teuchos::null)
768  Ap_ = MVT::Clone( *rhsMV, 1 );
769 
770  // Generate r_ only if it doesn't exist
771  if (r_ == Teuchos::null)
772  r_ = MVT::Clone( *rhsMV, 1 );
773 
774  // Generate z_ only if it doesn't exist
775  if (z_ == Teuchos::null)
776  z_ = MVT::Clone( *rhsMV, 1 );
777 
778  // If the recycle space has not been initialized before, generate it using the RHS from lp_.
779  if (U_ == Teuchos::null) {
780  U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
781  }
782  else {
783  // Generate U_ by cloning itself ONLY if more space is needed.
784  if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
785  Teuchos::RCP<const MV> tmp = U_;
786  U_ = MVT::Clone( *tmp, recycleBlocks_ );
787  }
788  }
789 
790  // If the recycle space has not be initialized before, generate it using the RHS from lp_.
791  if (AU_ == Teuchos::null) {
792  AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
793  }
794  else {
795  // Generate AU_ by cloning itself ONLY if more space is needed.
796  if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
797  Teuchos::RCP<const MV> tmp = AU_;
798  AU_ = MVT::Clone( *tmp, recycleBlocks_ );
799  }
800  }
801 
802  // If the recycle space has not been initialized before, generate it using the RHS from lp_.
803  if (U1_ == Teuchos::null) {
804  U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
805  }
806  else {
807  // Generate U1_ by cloning itself ONLY if more space is needed.
808  if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
809  Teuchos::RCP<const MV> tmp = U1_;
810  U1_ = MVT::Clone( *tmp, recycleBlocks_ );
811  }
812  }
813 
814  // Generate Alpha_ only if it doesn't exist, otherwise resize it.
815  if (Alpha_ == Teuchos::null)
816  Alpha_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
817  else {
818  if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
819  Alpha_->reshape( numBlocks_, 1 );
820  }
821 
822  // Generate Beta_ only if it doesn't exist, otherwise resize it.
823  if (Beta_ == Teuchos::null)
824  Beta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
825  else {
826  if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
827  Beta_->reshape( numBlocks_ + 1, 1 );
828  }
829 
830  // Generate D_ only if it doesn't exist, otherwise resize it.
831  if (D_ == Teuchos::null)
832  D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
833  else {
834  if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
835  D_->reshape( numBlocks_, 1 );
836  }
837 
838  // Generate Delta_ only if it doesn't exist, otherwise resize it.
839  if (Delta_ == Teuchos::null)
840  Delta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
841  else {
842  if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
843  Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
844  }
845 
846  // Generate UTAU_ only if it doesn't exist, otherwise resize it.
847  if (UTAU_ == Teuchos::null)
848  UTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
849  else {
850  if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
851  UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
852  }
853 
854  // Generate LUUTAU_ only if it doesn't exist, otherwise resize it.
855  if (LUUTAU_ == Teuchos::null)
856  LUUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
857  else {
858  if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
859  LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
860  }
861 
862  // Generate ipiv_ only if it doesn't exist, otherwise resize it.
863  if (ipiv_ == Teuchos::null)
864  ipiv_ = Teuchos::rcp( new std::vector<int>(recycleBlocks_) );
865  else {
866  if ( (int)ipiv_->size() != recycleBlocks_ ) // if ipiv not correct size, always resize it
867  ipiv_->resize(recycleBlocks_);
868  }
869 
870  // Generate AUTAU_ only if it doesn't exist, otherwise resize it.
871  if (AUTAU_ == Teuchos::null)
872  AUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
873  else {
874  if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
875  AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
876  }
877 
878  // Generate rTz_old_ only if it doesn't exist
879  if (rTz_old_ == Teuchos::null)
881  else {
882  if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
883  rTz_old_->reshape( 1, 1 );
884  }
885 
886  // Generate F_ only if it doesn't exist
887  if (F_ == Teuchos::null)
888  F_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
889  else {
890  if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
891  F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
892  }
893 
894  // Generate G_ only if it doesn't exist
895  if (G_ == Teuchos::null)
896  G_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
897  else {
898  if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
899  G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
900  }
901 
902  // Generate Y_ only if it doesn't exist
903  if (Y_ == Teuchos::null)
904  Y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
905  else {
906  if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
907  Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
908  }
909 
910  // Generate L2_ only if it doesn't exist
911  if (L2_ == Teuchos::null)
912  L2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
913  else {
914  if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
915  L2_->reshape( numBlocks_+1, numBlocks_ );
916  }
917 
918  // Generate DeltaL2_ only if it doesn't exist
919  if (DeltaL2_ == Teuchos::null)
920  DeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
921  else {
922  if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
923  DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
924  }
925 
926  // Generate AU1TUDeltaL2_ only if it doesn't exist
927  if (AU1TUDeltaL2_ == Teuchos::null)
928  AU1TUDeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
929  else {
930  if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
931  AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
932  }
933 
934  // Generate AU1TAU1_ only if it doesn't exist
935  if (AU1TAU1_ == Teuchos::null)
936  AU1TAU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
937  else {
938  if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
939  AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
940  }
941 
942  // Generate GY_ only if it doesn't exist
943  if (GY_ == Teuchos::null)
944  GY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
945  else {
946  if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
947  GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
948  }
949 
950  // Generate AU1TU1_ only if it doesn't exist
951  if (AU1TU1_ == Teuchos::null)
952  AU1TU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
953  else {
954  if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
955  AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
956  }
957 
958  // Generate FY_ only if it doesn't exist
959  if (FY_ == Teuchos::null)
960  FY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
961  else {
962  if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
963  FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
964  }
965 
966  // Generate AU1TAP_ only if it doesn't exist
967  if (AU1TAP_ == Teuchos::null)
968  AU1TAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
969  else {
970  if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
971  AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
972  }
973 
974  // Generate APTAP_ only if it doesn't exist
975  if (APTAP_ == Teuchos::null)
976  APTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
977  else {
978  if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
979  APTAP_->reshape( numBlocks_, numBlocks_ );
980  }
981 
982  // If the subspace has not been initialized before, generate it using the RHS from lp_.
983  if (U1Y1_ == Teuchos::null) {
984  U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
985  }
986  else {
987  // Generate U1Y1_ by cloning itself ONLY if more space is needed.
988  if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
989  Teuchos::RCP<const MV> tmp = U1Y1_;
990  U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
991  }
992  }
993 
994  // If the subspace has not been initialized before, generate it using the RHS from lp_.
995  if (PY2_ == Teuchos::null) {
996  PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
997  }
998  else {
999  // Generate PY2_ by cloning itself ONLY if more space is needed.
1000  if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1001  Teuchos::RCP<const MV> tmp = PY2_;
1002  PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1003  }
1004  }
1005 
1006  // Generate AUTAP_ only if it doesn't exist
1007  if (AUTAP_ == Teuchos::null)
1008  AUTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
1009  else {
1010  if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
1011  AUTAP_->reshape( recycleBlocks_, numBlocks_ );
1012  }
1013 
1014  // Generate AU1TU_ only if it doesn't exist
1015  if (AU1TU_ == Teuchos::null)
1016  AU1TU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
1017  else {
1018  if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
1019  AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
1020  }
1021 
1022 
1023  }
1024 }
1025 
1026 template<class ScalarType, class MV, class OP>
1028 
1030  std::vector<int> index(recycleBlocks_);
1031  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1032  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1033 
1034  // Count of number of cycles performed on current rhs
1035  int cycle = 0;
1036 
1037  // Set the current parameters if they were not set before.
1038  // NOTE: This may occur if the user generated the solver manager with the default constructor and
1039  // then didn't set any parameters using setParameters().
1040  if (!params_Set_) {
1041  setParameters(Teuchos::parameterList(*getValidParameters()));
1042  }
1043 
1045  "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1047  "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1048  TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1050  "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1051 
1052  // Grab the preconditioning object
1053  Teuchos::RCP<OP> precObj;
1054  if (problem_->getLeftPrec() != Teuchos::null) {
1055  precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1056  }
1057  else if (problem_->getRightPrec() != Teuchos::null) {
1058  precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1059  }
1060 
1061  // Create indices for the linear systems to be solved.
1062  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1063  std::vector<int> currIdx(1);
1064  currIdx[0] = 0;
1065 
1066  // Inform the linear problem of the current linear system to solve.
1067  problem_->setLSIndex( currIdx );
1068 
1069  // Check the number of blocks and change them if necessary.
1070  ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1071  if (numBlocks_ > dim) {
1072  numBlocks_ = Teuchos::asSafe<int>(dim);
1073  params_->set("Num Blocks", numBlocks_);
1074  printer_->stream(Warnings) <<
1075  "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1076  " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1077  }
1078 
1079  // Initialize storage for all state variables
1080  initializeStateStorage();
1081 
1082  // Parameter list
1083  Teuchos::ParameterList plist;
1084  plist.set("Num Blocks",numBlocks_);
1085  plist.set("Recycled Blocks",recycleBlocks_);
1086 
1087  // Reset the status test.
1088  outputTest_->reset();
1089 
1090  // Assume convergence is achieved, then let any failed convergence set this to false.
1091  bool isConverged = true;
1092 
1093  // Compute AU = A*U, UTAU = U'*AU, AUTAU = (AU)'*(AU)
1094  if (existU_) {
1095  index.resize(recycleBlocks_);
1096  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1097  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1098  index.resize(recycleBlocks_);
1099  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1100  Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1101  // Initialize AU
1102  problem_->applyOp( *Utmp, *AUtmp );
1103  // Initialize UTAU
1104  Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1105  MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1106  // Initialize AUTAU ( AUTAU = AU'*(M\AU) )
1107  Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1108  if ( precObj != Teuchos::null ) {
1109  index.resize(recycleBlocks_);
1110  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1111  index.resize(recycleBlocks_);
1112  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1113  Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage
1114  OPT::Apply( *precObj, *AUtmp, *PCAU );
1115  MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1116  } else {
1117  MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1118  }
1119  }
1120 
1122  // RCG solver
1123 
1125  rcg_iter = Teuchos::rcp( new RCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
1126 
1127  // Enter solve() iterations
1128  {
1129 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1130  Teuchos::TimeMonitor slvtimer(*timerSolve_);
1131 #endif
1132 
1133  while ( numRHS2Solve > 0 ) {
1134 
1135  // Debugging output to tell use if recycle space exists and will be used
1136  if (printer_->isVerbosity( Debug ) ) {
1137  if (existU_) printer_->print( Debug, "Using recycle space generated from previous call to solve()." );
1138  else printer_->print( Debug, "No recycle space exists." );
1139  }
1140 
1141  // Reset the number of iterations.
1142  rcg_iter->resetNumIters();
1143 
1144  // Set the current number of recycle blocks and subspace dimension with the RCG iteration.
1145  rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1146 
1147  // Reset the number of calls that the status test output knows about.
1148  outputTest_->resetNumCalls();
1149 
1150  // indicate that updated recycle space has not yet been generated for this linear system
1151  existU1_ = false;
1152 
1153  // reset cycle count
1154  cycle = 0;
1155 
1156  // Get the current residual
1157  problem_->computeCurrResVec( &*r_ );
1158 
1159  // If U exists, find best soln over this space first
1160  if (existU_) {
1161  // Solve linear system UTAU * y = (U'*r)
1162  Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
1163  index.resize(recycleBlocks_);
1164  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1165  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1166  MVT::MvTransMv( one, *Utmp, *r_, Utr );
1167  Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1168  Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1169  LUUTAUtmp.assign(UTAUtmp);
1170  int info = 0;
1171  lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
1172  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1173  "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1174 
1175  // Update solution (x = x + U*y)
1176  MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1177 
1178  // Update residual ( r = r - AU*y )
1179  index.resize(recycleBlocks_);
1180  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1181  Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1182  MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1183  }
1184 
1185  if ( precObj != Teuchos::null ) {
1186  OPT::Apply( *precObj, *r_, *z_ );
1187  } else {
1188  z_ = r_;
1189  }
1190 
1191  // rTz_old = r'*z
1192  MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1193 
1194  if ( existU_ ) {
1195  // mu = UTAU\(AU'*z);
1196  Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
1197  index.resize(recycleBlocks_);
1198  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1199  Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1200  MVT::MvTransMv( one, *AUtmp, *z_, mu );
1201  Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1202  char TRANS = 'N';
1203  int info;
1204  lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
1205  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1206  "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1207  // p = z - U*mu;
1208  index.resize( 1 );
1209  index[0] = 0;
1210  Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1211  MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1212  MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1213  } else {
1214  // p = z;
1215  index.resize( 1 );
1216  index[0] = 0;
1217  Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1218  MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1219  }
1220 
1221  // Set the new state and initialize the solver.
1222  RCGIterState<ScalarType,MV> newstate;
1223 
1224  // Create RCP views here
1225  index.resize( numBlocks_+1 );
1226  for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1227  newstate.P = MVT::CloneViewNonConst( *P_, index );
1228  index.resize( recycleBlocks_ );
1229  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1230  newstate.U = MVT::CloneViewNonConst( *U_, index );
1231  index.resize( recycleBlocks_ );
1232  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1233  newstate.AU = MVT::CloneViewNonConst( *AU_, index );
1234  newstate.Alpha = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) );
1235  newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) );
1236  newstate.D = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) );
1237  newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1238  newstate.LUUTAU = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) );
1239  // assign the rest of the values in the struct
1240  newstate.curDim = 1; // We have initialized the first search vector
1241  newstate.Ap = Ap_;
1242  newstate.r = r_;
1243  newstate.z = z_;
1244  newstate.existU = existU_;
1245  newstate.ipiv = ipiv_;
1246  newstate.rTz_old = rTz_old_;
1247 
1248  rcg_iter->initialize(newstate);
1249 
1250  while(1) {
1251 
1252  // tell rcg_iter to iterate
1253  try {
1254  rcg_iter->iterate();
1255 
1257  //
1258  // check convergence first
1259  //
1261  if ( convTest_->getStatus() == Passed ) {
1262  // We have convergence
1263  break; // break from while(1){rcg_iter->iterate()}
1264  }
1266  //
1267  // check for maximum iterations
1268  //
1270  else if ( maxIterTest_->getStatus() == Passed ) {
1271  // we don't have convergence
1272  isConverged = false;
1273  break; // break from while(1){rcg_iter->iterate()}
1274  }
1276  //
1277  // check if cycle complete; update for next cycle
1278  //
1280  else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1281  // index into P_ of last search vector generated this cycle
1282  int lastp = -1;
1283  // index into Beta_ of last entry generated this cycle
1284  int lastBeta = -1;
1285  if (recycleBlocks_ > 0) {
1286  if (!existU_) {
1287  if (cycle == 0) { // No U, no U1
1288 
1289  Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
1290  Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
1291  Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1292  Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1293  Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1294  Ftmp.putScalar(zero);
1295  Gtmp.putScalar(zero);
1296  for (int ii=0;ii<numBlocks_;ii++) {
1297  Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1298  if (ii > 0) {
1299  Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1300  Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1301  }
1302  Ftmp(ii,ii) = Dtmp(ii,0);
1303  }
1304 
1305  // compute harmonic Ritz vectors
1306  Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
1307  getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1308 
1309  // U1 = [P(:,1:end-1)*Y];
1310  index.resize( numBlocks_ );
1311  for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1312  Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1313  index.resize( recycleBlocks_ );
1314  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1315  Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1316  MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1317 
1318  // Precompute some variables for next cycle
1319 
1320  // AU1TAU1 = Y'*G*Y;
1321  Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
1322  GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1323  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1324  AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1325 
1326  // AU1TU1 = Y'*F*Y;
1327  Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
1328  FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1329  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1330  AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1331 
1332  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
1333  // Must reinitialize AU1TAP; can become dense later
1334  AU1TAPtmp.putScalar(zero);
1335  // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
1336  ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1337  for (int ii=0; ii<recycleBlocks_; ++ii) {
1338  AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1339  }
1340 
1341  // indicate that updated recycle space now defined
1342  existU1_ = true;
1343 
1344  // Indicate the size of the P, Beta structures generated this cycle
1345  lastp = numBlocks_;
1346  lastBeta = numBlocks_-1;
1347 
1348  } // if (cycle == 0)
1349  else { // No U, but U1 guaranteed to exist now
1350 
1351  // Finish computation of subblocks
1352  // AU1TAP = AU1TAP * D(1);
1353  Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1354  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1355  AU1TAPtmp.scale(Dtmp(0,0));
1356 
1357  Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1358  Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1359  Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1360  APTAPtmp.putScalar(zero);
1361  for (int ii=0; ii<numBlocks_; ii++) {
1362  APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1363  if (ii > 0) {
1364  APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1365  APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1366  }
1367  }
1368 
1369  // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
1370  Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1371  Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1372  Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1373  Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1374  Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1375  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1376  F11.assign(AU1TU1tmp);
1377  F12.putScalar(zero);
1378  F21.putScalar(zero);
1379  F22.putScalar(zero);
1380  for(int ii=0;ii<numBlocks_;ii++) {
1381  F22(ii,ii) = Dtmp(ii,0);
1382  }
1383 
1384  // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
1385  Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1386  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1387  Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1388  Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1389  Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1390  Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1391  G11.assign(AU1TAU1tmp);
1392  G12.assign(AU1TAPtmp);
1393  // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1394  for (int ii=0;ii<recycleBlocks_;++ii)
1395  for (int jj=0;jj<numBlocks_;++jj)
1396  G21(jj,ii) = G12(ii,jj);
1397  G22.assign(APTAPtmp);
1398 
1399  // compute harmonic Ritz vectors
1400  Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1401  getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1402 
1403  // U1 = [U1 P(:,2:end-1)]*Y;
1404  index.resize( numBlocks_ );
1405  for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1406  Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1407  index.resize( recycleBlocks_ );
1408  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1409  Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1410  Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1411  index.resize( recycleBlocks_ );
1412  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1413  Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1414  index.resize( recycleBlocks_ );
1415  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1416  Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1417  Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1418  MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1419  MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1420  MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1421 
1422  // Precompute some variables for next cycle
1423 
1424  // AU1TAU1 = Y'*G*Y;
1425  Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1426  GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1427  AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1428 
1429  // AU1TAP = zeros(k,m);
1430  // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
1431  AU1TAPtmp.putScalar(zero);
1432  ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1433  for (int ii=0; ii<recycleBlocks_; ++ii) {
1434  AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1435  }
1436 
1437  // AU1TU1 = Y'*F*Y;
1438  Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1439  FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1440  //Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1441  AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1442 
1443  // Indicate the size of the P, Beta structures generated this cycle
1444  lastp = numBlocks_+1;
1445  lastBeta = numBlocks_;
1446 
1447  } // if (cycle != 1)
1448  } // if (!existU_)
1449  else { // U exists
1450  if (cycle == 0) { // No U1, but U exists
1451  Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1452  Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1453  Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1454  Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1455  APTAPtmp.putScalar(zero);
1456  for (int ii=0; ii<numBlocks_; ii++) {
1457  APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1458  if (ii > 0) {
1459  APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1460  APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1461  }
1462  }
1463 
1464  Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1465  L2tmp.putScalar(zero);
1466  for(int ii=0;ii<numBlocks_;ii++) {
1467  L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1468  L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1469  }
1470 
1471  // AUTAP = UTAU*Delta*L2;
1472  Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
1473  Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1474  Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1475  Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1476  DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1477  AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
1478 
1479  // F = [UTAU zeros(k,m); zeros(m,k) diag(D)];
1480  Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1481  Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1482  Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1483  Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1484  Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1485  F11.assign(UTAUtmp);
1486  F12.putScalar(zero);
1487  F21.putScalar(zero);
1488  F22.putScalar(zero);
1489  for(int ii=0;ii<numBlocks_;ii++) {
1490  F22(ii,ii) = Dtmp(ii,0);
1491  }
1492 
1493  // G = [AUTAU AUTAP; AUTAP' APTAP];
1494  Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1495  Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1496  Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1497  Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1498  Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1499  Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1500  G11.assign(AUTAUtmp);
1501  G12.assign(AUTAPtmp);
1502  // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1503  for (int ii=0;ii<recycleBlocks_;++ii)
1504  for (int jj=0;jj<numBlocks_;++jj)
1505  G21(jj,ii) = G12(ii,jj);
1506  G22.assign(APTAPtmp);
1507 
1508  // compute harmonic Ritz vectors
1509  Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1510  getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1511 
1512  // U1 = [U P(:,1:end-1)]*Y;
1513  index.resize( recycleBlocks_ );
1514  for (int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1515  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1516  index.resize( numBlocks_ );
1517  for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1518  Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1519  index.resize( recycleBlocks_ );
1520  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1521  Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1522  Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1523  index.resize( recycleBlocks_ );
1524  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1525  Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1526  Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1527  index.resize( recycleBlocks_ );
1528  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1529  Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1530  MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1531  MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1532  MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1533 
1534  // Precompute some variables for next cycle
1535 
1536  // AU1TAU1 = Y'*G*Y;
1537  Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1538  GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1539  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1540  AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1541 
1542  // AU1TU1 = Y'*F*Y;
1543  Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1544  FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1545  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1546  AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1547 
1548  // AU1TU = UTAU;
1549  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1550  AU1TUtmp.assign(UTAUtmp);
1551 
1552  // dold = D(end);
1553  dold = Dtmp(numBlocks_-1,0);
1554 
1555  // indicate that updated recycle space now defined
1556  existU1_ = true;
1557 
1558  // Indicate the size of the P, Beta structures generated this cycle
1559  lastp = numBlocks_;
1560  lastBeta = numBlocks_-1;
1561  }
1562  else { // Have U and U1
1563  Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1564  Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1565  Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1566  Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1567  for (int ii=0; ii<numBlocks_; ii++) {
1568  APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1569  if (ii > 0) {
1570  APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1571  APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1572  }
1573  }
1574 
1575  Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1576  for(int ii=0;ii<numBlocks_;ii++) {
1577  L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1578  L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1579  }
1580 
1581  // M(end,1) = dold*(-Beta(1)/Alpha(1));
1582  // AU1TAP = Y'*[AU1TU*Delta*L2; M];
1583  Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1584  Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1585  DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1586  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
1587  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1588  AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
1589  Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1590  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1591  AU1TAPtmp.putScalar(zero);
1592  AU1TAPtmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
1593  Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1594  ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1595  for(int ii=0;ii<recycleBlocks_;ii++) {
1596  AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1597  }
1598 
1599  // AU1TU = Y1'*AU1TU
1600  Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
1601  Y1TAU1TU.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
1602  AU1TUtmp.assign(Y1TAU1TU);
1603 
1604  // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
1605  Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1606  Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1607  Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1608  Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1609  Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1610  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1611  F11.assign(AU1TU1tmp);
1612  for(int ii=0;ii<numBlocks_;ii++) {
1613  F22(ii,ii) = Dtmp(ii,0);
1614  }
1615 
1616  // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
1617  Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1618  Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1619  Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1620  Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1621  Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1622  Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1623  G11.assign(AU1TAU1tmp);
1624  G12.assign(AU1TAPtmp);
1625  // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1626  for (int ii=0;ii<recycleBlocks_;++ii)
1627  for (int jj=0;jj<numBlocks_;++jj)
1628  G21(jj,ii) = G12(ii,jj);
1629  G22.assign(APTAPtmp);
1630 
1631  // compute harmonic Ritz vectors
1632  Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1633  getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1634 
1635  // U1 = [U1 P(:,2:end-1)]*Y;
1636  index.resize( numBlocks_ );
1637  for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1638  Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1639  index.resize( recycleBlocks_ );
1640  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1641  Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1642  index.resize( recycleBlocks_ );
1643  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1644  Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1645  index.resize( recycleBlocks_ );
1646  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1647  Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1648  MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1649  MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1650  MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1651 
1652  // Precompute some variables for next cycle
1653 
1654  // AU1TAU1 = Y'*G*Y;
1655  Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1656  GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1657  AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1658 
1659  // AU1TU1 = Y'*F*Y;
1660  Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1661  FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1662  AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1663 
1664  // dold = D(end);
1665  dold = Dtmp(numBlocks_-1,0);
1666 
1667  // Indicate the size of the P, Beta structures generated this cycle
1668  lastp = numBlocks_+1;
1669  lastBeta = numBlocks_;
1670 
1671  }
1672  }
1673  } // if (recycleBlocks_ > 0)
1674 
1675  // Cleanup after end of cycle
1676 
1677  // P = P(:,end-1:end);
1678  index.resize( 2 );
1679  index[0] = lastp-1; index[1] = lastp;
1680  Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1681  index[0] = 0; index[1] = 1;
1682  MVT::SetBlock(*Ptmp2,index,*P_);
1683 
1684  // Beta = Beta(end);
1685  (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1686 
1687  // Delta = Delta(:,end);
1688  if (existU_) { // Delta only initialized if U exists
1689  Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
1690  Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1691  mu1.assign(mu2);
1692  }
1693 
1694  // Now reinitialize state variables for next cycle
1695  newstate.P = Teuchos::null;
1696  index.resize( numBlocks_+1 );
1697  for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1698  newstate.P = MVT::CloneViewNonConst( *P_, index );
1699 
1700  newstate.Beta = Teuchos::null;
1701  newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
1702 
1703  newstate.Delta = Teuchos::null;
1704  newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1705 
1706  newstate.curDim = 1; // We have initialized the first search vector
1707 
1708  // Pass to iteration object
1709  rcg_iter->initialize(newstate);
1710 
1711  // increment cycle count
1712  cycle = cycle + 1;
1713 
1714  }
1716  //
1717  // we returned from iterate(), but none of our status tests Passed.
1718  // something is wrong, and it is probably our fault.
1719  //
1721  else {
1722  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1723  "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1724  }
1725  }
1726  catch (const StatusTestNaNError& e) {
1727  // A NaN was detected in the solver. Set the solution to zero and return unconverged.
1728  achievedTol_ = MT::one();
1729  Teuchos::RCP<MV> X = problem_->getLHS();
1730  MVT::MvInit( *X, SCT::zero() );
1731  printer_->stream(Warnings) << "Belos::RCGSolMgr::solve(): Warning! NaN has been detected!"
1732  << std::endl;
1733  return Unconverged;
1734  }
1735  catch (const std::exception &e) {
1736  printer_->stream(Errors) << "Error! Caught std::exception in RCGIter::iterate() at iteration "
1737  << rcg_iter->getNumIters() << std::endl
1738  << e.what() << std::endl;
1739  throw;
1740  }
1741  }
1742 
1743  // Inform the linear problem that we are finished with this block linear system.
1744  problem_->setCurrLS();
1745 
1746  // Update indices for the linear systems to be solved.
1747  numRHS2Solve -= 1;
1748  if ( numRHS2Solve > 0 ) {
1749  currIdx[0]++;
1750  // Set the next indices.
1751  problem_->setLSIndex( currIdx );
1752  }
1753  else {
1754  currIdx.resize( numRHS2Solve );
1755  }
1756 
1757  // Update the recycle space for the next linear system
1758  if (existU1_) { // be sure updated recycle space was created
1759  // U = U1
1760  index.resize(recycleBlocks_);
1761  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1762  MVT::SetBlock(*U1_,index,*U_);
1763  // Set flag indicating recycle space is now defined
1764  existU_ = true;
1765  if (numRHS2Solve > 0) { // also update AU, UTAU, and AUTAU
1766  // Free pointers in newstate
1767  newstate.P = Teuchos::null;
1768  newstate.Ap = Teuchos::null;
1769  newstate.r = Teuchos::null;
1770  newstate.z = Teuchos::null;
1771  newstate.U = Teuchos::null;
1772  newstate.AU = Teuchos::null;
1773  newstate.Alpha = Teuchos::null;
1774  newstate.Beta = Teuchos::null;
1775  newstate.D = Teuchos::null;
1776  newstate.Delta = Teuchos::null;
1777  newstate.LUUTAU = Teuchos::null;
1778  newstate.ipiv = Teuchos::null;
1779  newstate.rTz_old = Teuchos::null;
1780 
1781  // Reinitialize AU, UTAU, AUTAU
1782  index.resize(recycleBlocks_);
1783  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1784  Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1785  index.resize(recycleBlocks_);
1786  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1787  Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1788  // Initialize AU
1789  problem_->applyOp( *Utmp, *AUtmp );
1790  // Initialize UTAU
1791  Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1792  MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1793  // Initialize AUTAU ( AUTAU = AU'*(M\AU) )
1794  Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1795  if ( precObj != Teuchos::null ) {
1796  index.resize(recycleBlocks_);
1797  for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1798  index.resize(recycleBlocks_);
1799  for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1800  Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage
1801  OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1802  MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1803  } else {
1804  MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1805  }
1806  } // if (numRHS2Solve > 0)
1807 
1808  } // if (existU1)
1809  }// while ( numRHS2Solve > 0 )
1810 
1811  }
1812 
1813  // print final summary
1814  sTest_->print( printer_->stream(FinalSummary) );
1815 
1816  // print timing information
1817 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1818  // Calling summarize() can be expensive, so don't call unless the
1819  // user wants to print out timing details. summarize() will do all
1820  // the work even if it's passed a "black hole" output stream.
1821  if (verbosity_ & TimingDetails)
1822  Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1823 #endif
1824 
1825  // get iteration information for this solve
1826  numIters_ = maxIterTest_->getNumIters();
1827 
1828  // Save the convergence test value ("achieved tolerance") for this solve.
1829  {
1830  using Teuchos::rcp_dynamic_cast;
1831  typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1832  // testValues is nonnull and not persistent.
1833  const std::vector<MagnitudeType>* pTestValues =
1834  rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1835 
1836  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1837  "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1838  "method returned NULL. Please report this bug to the Belos developers.");
1839 
1840  TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1841  "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1842  "method returned a vector of length zero. Please report this bug to the "
1843  "Belos developers.");
1844 
1845  // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1846  // achieved tolerances for all vectors in the current solve(), or
1847  // just for the vectors from the last deflation?
1848  achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1849  }
1850 
1851  if (!isConverged) {
1852  return Unconverged; // return from RCGSolMgr::solve()
1853  }
1854  return Converged; // return from RCGSolMgr::solve()
1855 }
1856 
1857 // Compute the harmonic eigenpairs of the projected, dense system.
1858 template<class ScalarType, class MV, class OP>
1862  // order of F,G
1863  int n = F.numCols();
1864 
1865  // The LAPACK interface
1867 
1868  // Magnitude of harmonic Ritz values
1869  std::vector<MagnitudeType> w(n);
1870 
1871  // Sorted order of harmonic Ritz values
1872  std::vector<int> iperm(n);
1873 
1874  // Compute k smallest harmonic Ritz pairs
1875  // SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO )
1876  int itype = 1; // solve A*x = (lambda)*B*x
1877  char jobz='V'; // compute eigenvalues and eigenvectors
1878  char uplo='U'; // since F,G symmetric, reference only their upper triangular data
1879  std::vector<ScalarType> work(1);
1880  int lwork = -1;
1881  int info = 0;
1882  // since SYGV destroys workspace, create copies of F,G
1885 
1886  // query for optimal workspace size
1887  lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1889  "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1890  lwork = (int)work[0];
1891  work.resize(lwork);
1892  lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1894  "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
1895 
1896 
1897  // Construct magnitude of each harmonic Ritz value
1898  this->sort(w,n,iperm);
1899 
1900  // Select recycledBlocks_ smallest eigenvectors
1901  for( int i=0; i<recycleBlocks_; i++ ) {
1902  for( int j=0; j<n; j++ ) {
1903  Y(j,i) = G2(j,iperm[i]);
1904  }
1905  }
1906 
1907 }
1908 
1909 // This method sorts list of n floating-point numbers and return permutation vector
1910 template<class ScalarType, class MV, class OP>
1911 void RCGSolMgr<ScalarType,MV,OP,true>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm)
1912 {
1913  int l, r, j, i, flag;
1914  int RR2;
1915  double dRR, dK;
1916 
1917  // Initialize the permutation vector.
1918  for(j=0;j<n;j++)
1919  iperm[j] = j;
1920 
1921  if (n <= 1) return;
1922 
1923  l = n / 2 + 1;
1924  r = n - 1;
1925  l = l - 1;
1926  dRR = dlist[l - 1];
1927  dK = dlist[l - 1];
1928 
1929  RR2 = iperm[l - 1];
1930  while (r != 0) {
1931  j = l;
1932  flag = 1;
1933 
1934  while (flag == 1) {
1935  i = j;
1936  j = j + j;
1937 
1938  if (j > r + 1)
1939  flag = 0;
1940  else {
1941  if (j < r + 1)
1942  if (dlist[j] > dlist[j - 1]) j = j + 1;
1943 
1944  if (dlist[j - 1] > dK) {
1945  dlist[i - 1] = dlist[j - 1];
1946  iperm[i - 1] = iperm[j - 1];
1947  }
1948  else {
1949  flag = 0;
1950  }
1951  }
1952  }
1953  dlist[i - 1] = dRR;
1954  iperm[i - 1] = RR2;
1955  if (l == 1) {
1956  dRR = dlist [r];
1957  RR2 = iperm[r];
1958  dK = dlist[r];
1959  dlist[r] = dlist[0];
1960  iperm[r] = iperm[0];
1961  r = r - 1;
1962  }
1963  else {
1964  l = l - 1;
1965  dRR = dlist[l - 1];
1966  RR2 = iperm[l - 1];
1967  dK = dlist[l - 1];
1968  }
1969  }
1970  dlist[0] = dRR;
1971  iperm[0] = RR2;
1972 }
1973 
1974 // This method requires the solver manager to return a std::string that describes itself.
1975 template<class ScalarType, class MV, class OP>
1977 {
1978  std::ostringstream oss;
1979  oss << "Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1980  return oss.str();
1981 }
1982 
1983 } // end Belos namespace
1984 
1985 #endif /* BELOS_RCG_SOLMGR_HPP */
ScalarType * values() const
Teuchos::RCP< MV > r
The current residual.
Collection of types and exceptions used within the Belos solvers.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Teuchos::RCP< MV > P
The current Krylov basis.
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
Teuchos::RCP< MV > AU
int getNumIters() const override
Get the iteration count for the most recent call to solve().
T & get(const std::string &name, T def_value)
static RCP< Time > getNewCounter(const std::string &name)
bool is_null(const std::shared_ptr< T > &p)
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Teuchos::RCP< MV > U
The recycled subspace and its image.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U.
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.
Implementation of the RCG (Recycling Conjugate Gradient) iterative linear solver. ...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
This class implements the RCG iteration, where a single-std::vector Krylov subspace is constructed...
An implementation of StatusTestResNorm using a family of residual norms.
int scale(const ScalarType alpha)
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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()
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
A factory class for generating StatusTestOutput objects.
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
bool isParameter(const std::string &name) const
Structure to contain pointers to RCGIter state variables.
Belos concrete class for performing the RCG iteration.
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:174
bool existU
Flag to indicate the recycle space should be used.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
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. ...
Teuchos::RCP< MV > z
The current preconditioned residual.
Teuchos::RCP< MV > Ap
A times current search vector.
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)
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
void SYGV(const OrdinalType &itype, const char &JOBZ, const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb, ScalarType *W, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
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::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
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 GESV(const OrdinalType &n, const OrdinalType &nrhs, ScalarType *A, const OrdinalType &lda, OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
void GETRS(const char &TRANS, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *A, const OrdinalType &lda, const OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
int curDim
The current dimension of the reduction.
RCGSolMgrLAPACKFailure(const std::string &what_arg)
bool isType(const std::string &name) const
RCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
A class for extending the status testing capabilities of Belos via logical combinations.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Class which defines basic traits for the operator type.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Belos header file which uses auto-configuration information to include necessary C++ headers...
RCGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
OrdinalType stride() const
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...

Generated on Fri Dec 6 2024 09:24:06 for Belos by doxygen 1.8.5