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"
24 #include "BelosStatusTestMaxIters.hpp"
25 #include "BelosStatusTestGenResNorm.hpp"
26 #include "BelosStatusTestCombo.hpp"
27 #include "BelosStatusTestOutputFactory.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 =
113  Belos::Details::LapackSupportsScalar<ScalarType>::value &&
114  ! Teuchos::ScalarTraits<ScalarType>::isComplex>
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 =
121  Belos::Details::LapackSupportsScalar<ScalarType>::value &&
122  ! Teuchos::ScalarTraits<ScalarType>::isComplex;
123  typedef Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP,
124  scalarTypeIsSupported> base_type;
125 
126  public:
127  RCGSolMgr () :
128  base_type ()
129  {}
130  RCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
131  const Teuchos::RCP<Teuchos::ParameterList> &pl) :
132  base_type ()
133  {}
134  virtual ~RCGSolMgr () {}
135 
137  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
138  return Teuchos::rcp(new RCGSolMgr<ScalarType,MV,OP,supportsScalarType>);
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:
149  typedef MultiVecTraits<ScalarType,MV> MVT;
150  typedef OperatorTraits<ScalarType,MV,OP> OPT;
151  typedef Teuchos::ScalarTraits<ScalarType> SCT;
152  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
153  typedef Teuchos::ScalarTraits<MagnitudeType> MT;
154 
155  public:
156 
158 
159 
165  RCGSolMgr();
166 
188  RCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
189  const Teuchos::RCP<Teuchos::ParameterList> &pl );
190 
192  virtual ~RCGSolMgr() {};
193 
195  Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
196  return Teuchos::rcp(new RCGSolMgr<ScalarType,MV,OP>);
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 
211  Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
212 
218  Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
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,
305  const Teuchos::SerialDenseMatrix<int,ScalarType> &G,
306  Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
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.
315  Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
316 
317  // Output manager.
318  Teuchos::RCP<OutputManager<ScalarType> > printer_;
319  Teuchos::RCP<std::ostream> outputStream_;
320 
321  // Status test.
322  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
323  Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
324  Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
325  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
326 
327  // Current parameter list.
328  Teuchos::RCP<Teuchos::ParameterList> params_;
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
392  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_;
393  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_;
394  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
395  //
396  // Solutions to local least-squares problems
397  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
398  //
399  // The matrix U^T A U
400  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_;
401  //
402  // LU factorization of U^T A U
403  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
404  //
405  // Data from LU factorization of UTAU
406  Teuchos::RCP<std::vector<int> > ipiv_;
407  //
408  // The matrix (AU)^T AU
409  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_;
410  //
411  // The scalar r'*z
412  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
413  //
414  // Matrices needed for calculation of harmonic Ritz eigenproblem
415  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_;
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_;
420  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_;
421  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_;
422  Teuchos::RCP<MV> U1Y1_, PY2_;
423  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_;
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>
438 RCGSolMgr<ScalarType,MV,OP,true>::RCGSolMgr():
439  achievedTol_(0.0),
440  numIters_(0)
441 {
442  init();
443 }
444 
445 // Basic Constructor
446 template<class ScalarType, class MV, class OP>
447 RCGSolMgr<ScalarType,MV,OP,true>::RCGSolMgr(
448  const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
449  const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
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>
465 void RCGSolMgr<ScalarType,MV,OP,true>::init()
466 {
467  outputStream_ = Teuchos::rcpFromRef(std::cout);
468  convtol_ = DefaultSolverParameters::convTol;
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>
516 void RCGSolMgr<ScalarType,MV,OP,true>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
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>
696 Teuchos::RCP<const Teuchos::ParameterList>
697 RCGSolMgr<ScalarType,MV,OP,true>::getValidParameters() const
698 {
699  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
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>
740 void RCGSolMgr<ScalarType,MV,OP,true>::initializeStateStorage() {
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)
880  rTz_old_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
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>
1027 ReturnType RCGSolMgr<ScalarType,MV,OP,true>::solve() {
1028 
1029  Teuchos::LAPACK<int,ScalarType> lapack;
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 
1044  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,RCGSolMgrLinearProblemFailure,
1045  "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1046  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),RCGSolMgrLinearProblemFailure,
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),
1049  RCGSolMgrLinearProblemFailure,
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 
1124  Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
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>
1859 void RCGSolMgr<ScalarType,MV,OP,true>::getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType>& F,
1860  const Teuchos::SerialDenseMatrix<int,ScalarType>& /* G */,
1861  Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
1862  // order of F,G
1863  int n = F.numCols();
1864 
1865  // The LAPACK interface
1866  Teuchos::LAPACK<int,ScalarType> lapack;
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
1883  Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_ );
1884  Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *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);
1888  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
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);
1893  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
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>
1976 std::string RCGSolMgr<ScalarType,MV,OP,true>::description() const
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 */
Teuchos::SerialDenseMatrix::values
ScalarType * values() const
Belos::RCGIterState::r
Teuchos::RCP< MV > r
The current residual.
Definition: BelosRCGIter.hpp:74
BelosTypes.hpp
Collection of types and exceptions used within the Belos solvers.
Belos::RCGIterState::rTz_old
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > rTz_old
Definition: BelosRCGIter.hpp:90
Belos::RCGSolMgr< ScalarType, MV, OP, true >::getCurrentParameters
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Definition: BelosRCGSolMgr.hpp:211
Belos::OutputManager
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Definition: BelosIteration.hpp:32
Belos::RCGIterState::P
Teuchos::RCP< MV > P
The current Krylov basis.
Definition: BelosRCGIter.hpp:68
Belos::RCGIterState::ipiv
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Definition: BelosRCGIter.hpp:98
BelosOutputManager.hpp
Class which manages the output and verbosity of the Belos solvers.
is_null
bool is_null(const boost::shared_ptr< T > &p)
Belos::RCGIterState::AU
Teuchos::RCP< MV > AU
Definition: BelosRCGIter.hpp:83
Belos::RCGSolMgr< ScalarType, MV, OP, true >::getNumIters
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Definition: BelosRCGSolMgr.hpp:231
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Teuchos::TimeMonitor::getNewCounter
static RCP< Time > getNewCounter(const std::string &name)
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
Teuchos::SerialDenseMatrix::multiply
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Belos::RCGIterState::U
Teuchos::RCP< MV > U
The recycled subspace and its image.
Definition: BelosRCGIter.hpp:83
Belos::RCGIterState::LUUTAU
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > LUUTAU
The LU factorization of the matrix U^T A U.
Definition: BelosRCGIter.hpp:96
Belos::Details::SolverManagerRequiresRealLapack
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
Definition: BelosSolverManager.hpp:409
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Belos::StatusTestOutputFactory
A factory class for generating StatusTestOutput objects.
Definition: BelosStatusTestOutputFactory.hpp:38
Belos::RCGSolMgr
Implementation of the RCG (Recycling Conjugate Gradient) iterative linear solver. ...
Definition: BelosRCGSolMgr.hpp:115
Belos::RCGIterState::Beta
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Beta
Definition: BelosRCGIter.hpp:88
Belos::RCGIter
This class implements the RCG iteration, where a single-std::vector Krylov subspace is constructed...
Definition: BelosRCGIter.hpp:113
Belos::StatusTestGenResNorm
An implementation of StatusTestResNorm using a family of residual norms.
Definition: BelosStatusTestGenResNorm.hpp:47
Teuchos::SerialDenseMatrix::scale
int scale(const ScalarType alpha)
Belos::RCGSolMgr< ScalarType, MV, OP, true >::clone
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Definition: BelosRCGSolMgr.hpp:195
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Belos::DefaultSolverParameters::convTol
static const double convTol
Default convergence tolerance.
Definition: BelosTypes.hpp:261
BelosStatusTestMaxIters.hpp
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::ScalarTraits::name
static std::string name()
Belos::RCGSolMgr< ScalarType, MV, OP, true >::isLOADetected
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Definition: BelosRCGSolMgr.hpp:236
BelosStatusTestOutputFactory.hpp
A factory class for generating StatusTestOutput objects.
Belos::RCGSolMgrLinearProblemFailure::RCGSolMgrLinearProblemFailure
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
Definition: BelosRCGSolMgr.hpp:93
Belos::RCGSolMgrLinearProblemFailure
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Definition: BelosRCGSolMgr.hpp:92
Belos::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: BelosMultiVecTraits.hpp:97
BelosStatusTestCombo.hpp
Belos::StatusTest for logically combining several status tests.
Teuchos::ParameterList::isParameter
bool isParameter(const std::string &name) const
Belos::RCGIterState
Structure to contain pointers to RCGIter state variables.
Definition: BelosRCGIter.hpp:60
BelosRCGIter.hpp
Belos concrete class for performing the RCG iteration.
Teuchos_LAPACK.hpp
Belos::StatusTestMaxIters
A Belos::StatusTest class for specifying a maximum number of iterations.
Definition: BelosStatusTestMaxIters.hpp:31
Belos::ResetType
ResetType
How to reset the solver.
Definition: BelosTypes.hpp:174
Belos::RCGIterState::existU
bool existU
Flag to indicate the recycle space should be used.
Definition: BelosRCGIter.hpp:80
Belos::RCGSolMgr< ScalarType, MV, OP, true >::getProblem
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Definition: BelosRCGSolMgr.hpp:203
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
BelosSolverManager.hpp
Pure virtual base class which describes the basic interface for a solver manager. ...
Belos::RCGIterState::z
Teuchos::RCP< MV > z
The current preconditioned residual.
Definition: BelosRCGIter.hpp:77
Belos::RCGIterState::Ap
Teuchos::RCP< MV > Ap
A times current search vector.
Definition: BelosRCGIter.hpp:71
Teuchos::TimeMonitor::summarize
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)
Teuchos::SerialDenseMatrix::putScalar
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Belos::LinearProblem
A linear system to solve, and its associated information.
Definition: BelosIteration.hpp:29
BelosLinearProblem.hpp
Class which describes the linear problem to be solved by the iterative solver.
Teuchos::LAPACK::SYGV
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
Belos::Details::LapackSupportsScalar
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
Definition: BelosSolverManager.hpp:262
Belos::ReturnType
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:123
Teuchos::ParameterList::validateParameters
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Belos::RCGSolMgr< ScalarType, MV, OP, true >::getTimers
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Definition: BelosRCGSolMgr.hpp:218
Belos::RCGIterState::Alpha
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Alpha
Coefficients arising in RCG iteration.
Definition: BelosRCGIter.hpp:87
Teuchos::SerialDenseMatrix::numCols
OrdinalType numCols() const
Belos::StatusTestOutputFactory::create
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.
Definition: BelosStatusTestOutputFactory.hpp:79
Teuchos::LAPACK::GESV
void GESV(const OrdinalType &n, const OrdinalType &nrhs, ScalarType *A, const OrdinalType &lda, OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
Teuchos::LAPACK::GETRS
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
BelosStatusTestGenResNorm.hpp
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::RCGIterState::curDim
int curDim
The current dimension of the reduction.
Definition: BelosRCGIter.hpp:65
Belos::RCGSolMgrLAPACKFailure::RCGSolMgrLAPACKFailure
RCGSolMgrLAPACKFailure(const std::string &what_arg)
Definition: BelosRCGSolMgr.hpp:103
Teuchos::ParameterList::isType
bool isType(const std::string &name) const
Belos::RCGSolMgr::RCGSolMgr
RCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Definition: BelosRCGSolMgr.hpp:130
Belos::StatusTestCombo
A class for extending the status testing capabilities of Belos via logical combinations.
Definition: BelosStatusTestCombo.hpp:59
Belos::RCGSolMgr< ScalarType, MV, OP, true >::achievedTol
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Definition: BelosRCGSolMgr.hpp:226
Belos::OperatorTraits
Class which defines basic traits for the operator type.
Definition: BelosOperatorTraits.hpp:77
Belos::RCGIterState::D
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > D
Definition: BelosRCGIter.hpp:89
Teuchos_TimeMonitor.hpp
Belos::BelosError
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:28
Belos::RCGSolMgr::clone
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Definition: BelosRCGSolMgr.hpp:137
BelosConfigDefs.hpp
Belos header file which uses auto-configuration information to include necessary C++ headers...
Belos::RCGSolMgrLAPACKFailure
RCGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
Definition: BelosRCGSolMgr.hpp:102
Teuchos::SerialDenseMatrix::assign
SerialDenseMatrix< OrdinalType, ScalarType > & assign(const SerialDenseMatrix< OrdinalType, ScalarType > &Source)
Teuchos_as.hpp
Belos::RCGIterState::Delta
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > Delta
Solutions to local least-squares problems.
Definition: BelosRCGIter.hpp:93
Belos::RCGSolMgr< ScalarType, MV, OP, true >::setProblem
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
Definition: BelosRCGSolMgr.hpp:244
Teuchos::SerialDenseMatrix::stride
OrdinalType stride() const
Belos::RCGSolMgr< ScalarType, MV, OP, true >::reset
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Definition: BelosRCGSolMgr.hpp:259
Generated on Wed Oct 23 2024 09:24:27 for Belos by doxygen 1.8.5