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