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