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

Generated on Fri Aug 14 2020 10:48:34 for Belos by doxygen 1.8.5