BelosRCGSolMgr.hpp

Go to the documentation of this file.

//@HEADER
// ************************************************************************
//
//                 Belos: Block Linear Solvers Package
//                  Copyright 2004 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER

#ifndef BELOS_RCG_SOLMGR_HPP
#define BELOS_RCG_SOLMGR_HPP

#include "BelosConfigDefs.hpp"
#include "BelosTypes.hpp"

#include "BelosLinearProblem.hpp"
#include "BelosSolverManager.hpp"

#include "BelosRCGIter.hpp"
#include "BelosDGKSOrthoManager.hpp"
#include "BelosICGSOrthoManager.hpp"
#include "BelosIMGSOrthoManager.hpp"
#include "BelosStatusTestMaxIters.hpp"
#include "BelosStatusTestGenResNorm.hpp"
#include "BelosStatusTestCombo.hpp"
#include "BelosStatusTestOutputFactory.hpp"
#include "BelosOutputManager.hpp"
#include "Teuchos_BLAS.hpp"
#include "Teuchos_LAPACK.hpp"
#ifdef BELOS_TEUCHOS_TIME_MONITOR
#include "Teuchos_TimeMonitor.hpp"
#endif

namespace Belos {
  


  class RCGSolMgrLinearProblemFailure : public BelosError {public:
    RCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
    {}};
 
  class RCGSolMgrLAPACKFailure : public BelosError {public:
    RCGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
    {}};

  class RCGSolMgrRecyclingFailure : public BelosError {public:
    RCGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg)
    {}};

  
  template<class ScalarType, class MV, class OP>
  class RCGSolMgr : public SolverManager<ScalarType,MV,OP> {
    
  private:
    typedef MultiVecTraits<ScalarType,MV> MVT;
    typedef OperatorTraits<ScalarType,MV,OP> OPT;
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
    typedef Teuchos::ScalarTraits<MagnitudeType> MT;
    
  public:
    


     RCGSolMgr();

    RCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
       const Teuchos::RCP<Teuchos::ParameterList> &pl );
    
    virtual ~RCGSolMgr() {};
    

    
    const LinearProblem<ScalarType,MV,OP>& getProblem() const {
      return *problem_;
    }

    Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
    
    Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
    
    Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
      return Teuchos::tuple(timerSolve_);
    }

    MagnitudeType achievedTol() const {
      return achievedTol_;
    }

    int getNumIters() const {
      return numIters_;
    }

    bool isLOADetected() const { return false; }
 
    

   
    void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
   
    void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
    
   


    void reset( const ResetType type ) { 
      if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); 
      else if (type & Belos::RecycleSubspace) existU_ = false;
    }
 

    
    ReturnType solve();
    
    
    
    std::string description() const;
    
    
  private:

    // Called by all constructors; Contains init instructions common to all constructors
    void init();
 
    //  Computes harmonic eigenpairs of projected matrix created during one cycle.
    //  Y contains the harmonic Ritz vectors corresponding to the recycleBlocks eigenvalues of smallest magnitude.
    void getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType> &F, 
                         const Teuchos::SerialDenseMatrix<int,ScalarType> &G, 
                         Teuchos::SerialDenseMatrix<int,ScalarType>& Y);
 
    // Sort list of n floating-point numbers and return permutation vector
    void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);

    // Initialize solver state storage
    void initializeStateStorage();

    // Linear problem.
    Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
    
    // Output manager.
    Teuchos::RCP<OutputManager<ScalarType> > printer_;
    Teuchos::RCP<std::ostream> outputStream_;

    // Status test.
    Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
    Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
    Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
    Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;

    // Current parameter list.
    Teuchos::RCP<Teuchos::ParameterList> params_;
    
    // Default solver values.
    static const MagnitudeType convtol_default_;
    static const int maxIters_default_;
    static const int blockSize_default_;
    static const int numBlocks_default_;
    static const int recycleBlocks_default_;
    static const bool showMaxResNormOnly_default_;
    static const int verbosity_default_;
    static const int outputStyle_default_;
    static const int outputFreq_default_;
    static const std::string label_default_;
    static const Teuchos::RCP<std::ostream> outputStream_default_;

    //
    // Current solver values.
    //

    MagnitudeType convtol_;

    MagnitudeType achievedTol_;

    int maxIters_;

    int numIters_;

    int numBlocks_, recycleBlocks_;
    bool showMaxResNormOnly_;
    int verbosity_, outputStyle_, outputFreq_;

    // Solver State Storage
    // Search vectors
    Teuchos::RCP<MV> P_;
    //
    // A times current search direction
    Teuchos::RCP<MV> Ap_;
    //
    // Residual vector
    Teuchos::RCP<MV> r_;
    //
    // Preconditioned residual
    Teuchos::RCP<MV> z_;
    //
    // Flag indicating that the recycle space should be used
    bool existU_;
    //
    // Flag indicating that the updated recycle space has been created
    bool existU1_;
    //
    // Recycled subspace and its image
    Teuchos::RCP<MV> U_, AU_;
    //
    // Recycled subspace for next system and its image
    Teuchos::RCP<MV> U1_;
    //
    // Coefficients arising in RCG iteration
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_;
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_;
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
    //
    // Solutions to local least-squares problems
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_;
    //
    // The matrix U^T A U
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_;
    //
    // LU factorization of U^T A U
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_;
    //
    // Data from LU factorization of UTAU
    Teuchos::RCP<std::vector<int> > ipiv_;
    //
    // The matrix (AU)^T AU
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_;
    //
    // The scalar r'*z
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_;
    //
    // Matrices needed for calculation of harmonic Ritz eigenproblem
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_;
    //
    // Matrices needed for updating recycle space
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_;
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_;
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_;
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_;
    Teuchos::RCP<MV> U1Y1_, PY2_;
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_;
    ScalarType dold;
    
    // Timers.
    std::string label_;
    Teuchos::RCP<Teuchos::Time> timerSolve_;

    // Internal state variables.
    bool params_Set_;
  };


// Default solver values.
template<class ScalarType, class MV, class OP>
const typename RCGSolMgr<ScalarType,MV,OP>::MagnitudeType RCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;

template<class ScalarType, class MV, class OP>
const int RCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;

template<class ScalarType, class MV, class OP>
const int RCGSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 25;

template<class ScalarType, class MV, class OP>
const int RCGSolMgr<ScalarType,MV,OP>::blockSize_default_ = 1;
 
template<class ScalarType, class MV, class OP>
const int RCGSolMgr<ScalarType,MV,OP>::recycleBlocks_default_ = 3;

template<class ScalarType, class MV, class OP>
const bool RCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;

template<class ScalarType, class MV, class OP>
const int RCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;

template<class ScalarType, class MV, class OP>
const int RCGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;

template<class ScalarType, class MV, class OP>
const int RCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;

template<class ScalarType, class MV, class OP>
const std::string RCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";

template<class ScalarType, class MV, class OP>
const Teuchos::RCP<std::ostream> RCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);

// Empty Constructor
template<class ScalarType, class MV, class OP>
RCGSolMgr<ScalarType,MV,OP>::RCGSolMgr() {
  init();
}

// Basic Constructor
template<class ScalarType, class MV, class OP>
RCGSolMgr<ScalarType,MV,OP>::RCGSolMgr( 
                 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
                 const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
  problem_(problem)
{
  init();
  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");

  // If the parameter list pointer is null, then set the current parameters to the default parameter list.
  if ( !is_null(pl) ) {
    setParameters( pl );  
  }
}

// Common instructions executed in all constructors
template<class ScalarType, class MV, class OP>
void RCGSolMgr<ScalarType,MV,OP>::init()
{
  outputStream_ = outputStream_default_;
  convtol_ = convtol_default_;
  maxIters_ = maxIters_default_;
  numBlocks_ = numBlocks_default_;
  recycleBlocks_ = recycleBlocks_default_;
  verbosity_ = verbosity_default_;
  outputStyle_= outputStyle_default_;
  outputFreq_= outputFreq_default_;
  showMaxResNormOnly_ = showMaxResNormOnly_default_;
  label_ = label_default_;
  params_Set_ = false;
  P_ = Teuchos::null;
  Ap_ = Teuchos::null;
  r_ = Teuchos::null;
  z_ = Teuchos::null;
  existU_ = false;
  existU1_ = false;
  U_ = Teuchos::null;
  AU_ = Teuchos::null;
  U1_ = Teuchos::null;
  Alpha_ = Teuchos::null;
  Beta_ = Teuchos::null;
  D_ = Teuchos::null;
  Delta_ = Teuchos::null;
  UTAU_ = Teuchos::null;
  LUUTAU_ = Teuchos::null;
  ipiv_ = Teuchos::null;
  AUTAU_ = Teuchos::null;
  rTz_old_ = Teuchos::null;
  F_ = Teuchos::null;
  G_ = Teuchos::null;
  Y_ = Teuchos::null;
  L2_ = Teuchos::null;
  DeltaL2_ = Teuchos::null;
  AU1TUDeltaL2_ = Teuchos::null;
  AU1TAU1_ = Teuchos::null;
  AU1TU1_ = Teuchos::null;
  AU1TAP_ = Teuchos::null;
  FY_ = Teuchos::null;
  GY_ = Teuchos::null;
  APTAP_ = Teuchos::null;
  U1Y1_ = Teuchos::null;
  PY2_ = Teuchos::null;
  AUTAP_ = Teuchos::null;
  AU1TU_ = Teuchos::null;
  dold = 0.;
}

template<class ScalarType, class MV, class OP>
void RCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
{
  // Create the internal parameter list if ones doesn't already exist.
  if (params_ == Teuchos::null) {
    params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
  }
  else {
    params->validateParameters(*getValidParameters());
  }

  // Check for maximum number of iterations
  if (params->isParameter("Maximum Iterations")) {
    maxIters_ = params->get("Maximum Iterations",maxIters_default_);

    // Update parameter in our list and in status test.
    params_->set("Maximum Iterations", maxIters_);
    if (maxIterTest_!=Teuchos::null)
      maxIterTest_->setMaxIters( maxIters_ );
  }

  // Check for the maximum number of blocks.
  if (params->isParameter("Num Blocks")) {
    numBlocks_ = params->get("Num Blocks",numBlocks_default_);
    TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
                       "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
 
    // Update parameter in our list.
    params_->set("Num Blocks", numBlocks_);
  }
 
  // Check for the maximum number of blocks.
  if (params->isParameter("Num Recycled Blocks")) {
    recycleBlocks_ = params->get("Num Recycled Blocks",recycleBlocks_default_);
    TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
                       "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
 
    TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
                       "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
 
    // Update parameter in our list.
    params_->set("Num Recycled Blocks", recycleBlocks_);
  }

  // Check to see if the timer label changed.
  if (params->isParameter("Timer Label")) {
    std::string tempLabel = params->get("Timer Label", label_default_);

    // Update parameter in our list and solver timer
    if (tempLabel != label_) {
      label_ = tempLabel;
      params_->set("Timer Label", label_);
      std::string solveLabel = label_ + ": RCGSolMgr total solve time";
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
#endif
    }
  }

  // Check for a change in verbosity level
  if (params->isParameter("Verbosity")) {
    if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
      verbosity_ = params->get("Verbosity", verbosity_default_);
    } else {
      verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
    }

    // Update parameter in our list.
    params_->set("Verbosity", verbosity_);
    if (printer_ != Teuchos::null)
      printer_->setVerbosity(verbosity_);
  }

  // Check for a change in output style
  if (params->isParameter("Output Style")) {
    if (Teuchos::isParameterType<int>(*params,"Output Style")) {
      outputStyle_ = params->get("Output Style", outputStyle_default_);
    } else {
      outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
    }

    // Reconstruct the convergence test if the explicit residual test is not being used.
    params_->set("Output Style", outputStyle_);
    outputTest_ = Teuchos::null;
  }

  // output stream
  if (params->isParameter("Output Stream")) {
    outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");

    // Update parameter in our list.
    params_->set("Output Stream", outputStream_);
    if (printer_ != Teuchos::null)
      printer_->setOStream( outputStream_ );
  }

  // frequency level
  if (verbosity_ & Belos::StatusTestDetails) {
    if (params->isParameter("Output Frequency")) {
      outputFreq_ = params->get("Output Frequency", outputFreq_default_);
    }

    // Update parameter in out list and output status test.
    params_->set("Output Frequency", outputFreq_);
    if (outputTest_ != Teuchos::null)
      outputTest_->setOutputFrequency( outputFreq_ );
  }

  // Create output manager if we need to.
  if (printer_ == Teuchos::null) {
    printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
  }  
  
  // Convergence
  typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
  typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;

  // Check for convergence tolerance
  if (params->isParameter("Convergence Tolerance")) {
    convtol_ = params->get("Convergence Tolerance",convtol_default_);

    // Update parameter in our list and residual tests.
    params_->set("Convergence Tolerance", convtol_);
    if (convTest_ != Teuchos::null)
      convTest_->setTolerance( convtol_ );
  }
  
  if (params->isParameter("Show Maximum Residual Norm Only")) {
    showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");

    // Update parameter in our list and residual tests
    params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
    if (convTest_ != Teuchos::null)
      convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
  }

  // Create status tests if we need to.

  // Basic test checks maximum iterations and native residual.
  if (maxIterTest_ == Teuchos::null)
    maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
  
  // Implicit residual test, using the native residual to determine if convergence was achieved.
  if (convTest_ == Teuchos::null)
    convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
  
  if (sTest_ == Teuchos::null)
    sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
  
  if (outputTest_ == Teuchos::null) {

    // Create the status test output class.
    // This class manages and formats the output from the status test.
    StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
    outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
    
    // Set the solver string for the output test
    std::string solverDesc = " Recycling CG ";
    outputTest_->setSolverDesc( solverDesc );
  }

  // Create the timer if we need to.
  if (timerSolve_ == Teuchos::null) {
    std::string solveLabel = label_ + ": RCGSolMgr total solve time";
#ifdef BELOS_TEUCHOS_TIME_MONITOR
    timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
#endif
  }

  // Inform the solver manager that the current parameters were set.
  params_Set_ = true;
}

    
template<class ScalarType, class MV, class OP>
Teuchos::RCP<const Teuchos::ParameterList> 
RCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
{
  static Teuchos::RCP<const Teuchos::ParameterList> validPL;
  
  // Set all the valid parameters and their default values.
  if(is_null(validPL)) {
    Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
    pl->set("Convergence Tolerance", convtol_default_,
      "The relative residual tolerance that needs to be achieved by the\n"
      "iterative solver in order for the linear system to be declared converged.");
    pl->set("Maximum Iterations", maxIters_default_,
      "The maximum number of block iterations allowed for each\n"
      "set of RHS solved.");
    pl->set("Block Size", blockSize_default_,
      "Block Size Parameter -- currently must be 1 for RCG");
    pl->set("Num Blocks", numBlocks_default_,
      "The length of a cycle (and this max number of search vectors kept)\n");
    pl->set("Num Recycled Blocks", recycleBlocks_default_,
      "The number of vectors in the recycle subspace.");
    pl->set("Verbosity", verbosity_default_,
      "What type(s) of solver information should be outputted\n"
      "to the output stream.");
    pl->set("Output Style", outputStyle_default_,
      "What style is used for the solver information outputted\n"
      "to the output stream.");
    pl->set("Output Frequency", outputFreq_default_,
      "How often convergence information should be outputted\n"
      "to the output stream.");  
    pl->set("Output Stream", outputStream_default_,
      "A reference-counted pointer to the output stream where all\n"
      "solver output is sent.");
    pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
      "When convergence information is printed, only show the maximum\n"
      "relative residual norm when the block size is greater than one.");
    pl->set("Timer Label", label_default_,
      "The string to use as a prefix for the timer labels.");
    //  pl->set("Restart Timers", restartTimers_);
    validPL = pl;
  }
  return validPL;
}

// initializeStateStorage
template<class ScalarType, class MV, class OP>
void RCGSolMgr<ScalarType,MV,OP>::initializeStateStorage() {

    // Check if there is any multivector to clone from.
    Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
    if (rhsMV == Teuchos::null) {
      // Nothing to do
      return;
    }
    else {

      // Initialize the state storage
      TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument,
                         "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");

      // If the subspace has not been initialized before, generate it using the RHS from lp_.
      if (P_ == Teuchos::null) {
        P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
      }
      else {
        // Generate P_ by cloning itself ONLY if more space is needed.
        if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
          Teuchos::RCP<const MV> tmp = P_;
          P_ = MVT::Clone( *tmp, numBlocks_+2 );
        }
      }

      // Generate Ap_ only if it doesn't exist
      if (Ap_ == Teuchos::null)
        Ap_ = MVT::Clone( *rhsMV, 1 );

      // Generate r_ only if it doesn't exist
      if (r_ == Teuchos::null)
        r_ = MVT::Clone( *rhsMV, 1 );

      // Generate z_ only if it doesn't exist
      if (z_ == Teuchos::null)
        z_ = MVT::Clone( *rhsMV, 1 );

      // If the recycle space has not been initialized before, generate it using the RHS from lp_.
      if (U_ == Teuchos::null) {
        U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
      }
      else {
        // Generate U_ by cloning itself ONLY if more space is needed.
        if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
          Teuchos::RCP<const MV> tmp = U_;
          U_ = MVT::Clone( *tmp, recycleBlocks_ );
        }
      }

      // If the recycle space has not be initialized before, generate it using the RHS from lp_.
      if (AU_ == Teuchos::null) {
        AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
      }
      else {
        // Generate AU_ by cloning itself ONLY if more space is needed.
        if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
          Teuchos::RCP<const MV> tmp = AU_;
          AU_ = MVT::Clone( *tmp, recycleBlocks_ );
        }
      }

      // If the recycle space has not been initialized before, generate it using the RHS from lp_.
      if (U1_ == Teuchos::null) {
        U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
      }
      else {
        // Generate U1_ by cloning itself ONLY if more space is needed.
        if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
          Teuchos::RCP<const MV> tmp = U1_;
          U1_ = MVT::Clone( *tmp, recycleBlocks_ );
        }
      }

      // Generate Alpha_ only if it doesn't exist, otherwise resize it.
      if (Alpha_ == Teuchos::null)
        Alpha_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) );
      else {
        if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) )
          Alpha_->reshape( numBlocks_, 1 );
      }

      // Generate Beta_ only if it doesn't exist, otherwise resize it.
      if (Beta_ == Teuchos::null)
        Beta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) );
      else {
        if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) )
          Beta_->reshape( numBlocks_ + 1, 1 );
      }

      // Generate D_ only if it doesn't exist, otherwise resize it.
      if (D_ == Teuchos::null)
        D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) );
      else {
        if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) )
          D_->reshape( numBlocks_, 1 );
      }

      // Generate Delta_ only if it doesn't exist, otherwise resize it.
      if (Delta_ == Teuchos::null)
        Delta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) );
      else {
        if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) )
          Delta_->reshape( recycleBlocks_, numBlocks_ + 1 );
      }

      // Generate UTAU_ only if it doesn't exist, otherwise resize it.
      if (UTAU_ == Teuchos::null)
        UTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) )
          UTAU_->reshape( recycleBlocks_, recycleBlocks_ );
      }

      // Generate LUUTAU_ only if it doesn't exist, otherwise resize it.
      if (LUUTAU_ == Teuchos::null)
        LUUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) )
          LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
      }

      // Generate ipiv_ only if it doesn't exist, otherwise resize it.
      if (ipiv_ == Teuchos::null)
        ipiv_ = Teuchos::rcp( new std::vector<int>(recycleBlocks_) );
      else {
        if ( (int)ipiv_->size() != recycleBlocks_ ) // if ipiv not correct size, always resize it
          ipiv_->resize(recycleBlocks_);
      }

      // Generate AUTAU_ only if it doesn't exist, otherwise resize it.
      if (AUTAU_ == Teuchos::null)
        AUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) )
          AUTAU_->reshape( recycleBlocks_, recycleBlocks_ );
      }

      // Generate rTz_old_ only if it doesn't exist
      if (rTz_old_ == Teuchos::null)
        rTz_old_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) );
      else {
        if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) )
          rTz_old_->reshape( 1, 1 );
      }

      // Generate F_ only if it doesn't exist
      if (F_ == Teuchos::null)
        F_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
      else {
        if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) )
          F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
      }

      // Generate G_ only if it doesn't exist
      if (G_ == Teuchos::null)
        G_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) );
      else {
        if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) )
          G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
      }

      // Generate Y_ only if it doesn't exist
      if (Y_ == Teuchos::null)
        Y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) )
          Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
      }

      // Generate L2_ only if it doesn't exist
      if (L2_ == Teuchos::null)
        L2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) );
      else {
        if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) )
          L2_->reshape( numBlocks_+1, numBlocks_ );
      }

      // Generate DeltaL2_ only if it doesn't exist
      if (DeltaL2_ == Teuchos::null)
        DeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
      else {
        if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) )
          DeltaL2_->reshape( recycleBlocks_, numBlocks_ );
      }

      // Generate AU1TUDeltaL2_ only if it doesn't exist
      if (AU1TUDeltaL2_ == Teuchos::null)
        AU1TUDeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
      else {
        if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) )
          AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ );
      }

      // Generate AU1TAU1_ only if it doesn't exist
      if (AU1TAU1_ == Teuchos::null)
        AU1TAU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) )
          AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ );
      }

      // Generate GY_ only if it doesn't exist
      if (GY_ == Teuchos::null)
        GY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) )
          GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
      }

      // Generate AU1TU1_ only if it doesn't exist
      if (AU1TU1_ == Teuchos::null)
        AU1TU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) )
          AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ );
      }

      // Generate FY_ only if it doesn't exist
      if (FY_ == Teuchos::null)
        FY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) )
          FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ );
      }

      // Generate AU1TAP_ only if it doesn't exist
      if (AU1TAP_ == Teuchos::null)
        AU1TAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
      else {
        if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) )
          AU1TAP_->reshape( recycleBlocks_, numBlocks_ );
      }

      // Generate APTAP_ only if it doesn't exist
      if (APTAP_ == Teuchos::null)
        APTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) );
      else {
        if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) )
          APTAP_->reshape( numBlocks_, numBlocks_ );
      }

      // If the subspace has not been initialized before, generate it using the RHS from lp_.
      if (U1Y1_ == Teuchos::null) {
        U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
      }
      else {
        // Generate U1Y1_ by cloning itself ONLY if more space is needed.
        if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
          Teuchos::RCP<const MV> tmp = U1Y1_;
          U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
        }
      }

      // If the subspace has not been initialized before, generate it using the RHS from lp_.
      if (PY2_ == Teuchos::null) {
        PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
      }
      else {
        // Generate PY2_ by cloning itself ONLY if more space is needed.
        if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
          Teuchos::RCP<const MV> tmp = PY2_;
          PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
        }
      }

      // Generate AUTAP_ only if it doesn't exist
      if (AUTAP_ == Teuchos::null)
        AUTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) );
      else {
        if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) )
          AUTAP_->reshape( recycleBlocks_, numBlocks_ );
      }

      // Generate AU1TU_ only if it doesn't exist
      if (AU1TU_ == Teuchos::null)
        AU1TU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) );
      else {
        if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) )
          AU1TU_->reshape( recycleBlocks_, recycleBlocks_ );
      }


    }
}

template<class ScalarType, class MV, class OP>
ReturnType RCGSolMgr<ScalarType,MV,OP>::solve() {

  Teuchos::LAPACK<int,ScalarType> lapack;
  std::vector<int> index(recycleBlocks_);
  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();

  // Count of number of cycles performed on current rhs
  int cycle = 0;

  // Set the current parameters if they were not set before.
  // NOTE:  This may occur if the user generated the solver manager with the default constructor and 
  // then didn't set any parameters using setParameters().
  if (!params_Set_) {
    setParameters(Teuchos::parameterList(*getValidParameters()));
  }

  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,RCGSolMgrLinearProblemFailure,
                     "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),RCGSolMgrLinearProblemFailure,
                     "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");

  // Create indices for the linear systems to be solved.
  int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
  std::vector<int> currIdx(1);
  currIdx[0] = 0;

  // Inform the linear problem of the current linear system to solve.
  problem_->setLSIndex( currIdx );

  // Check the number of blocks and change them if necessary.
  int dim = MVT::GetVecLength( *(problem_->getRHS()) );
  if (numBlocks_ > dim) {
    numBlocks_ = dim;
    params_->set("Num Blocks", numBlocks_);
    printer_->stream(Warnings) <<
      "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
      " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
  }

  // Initialize storage for all state variables
  initializeStateStorage();

  // Parameter list
  Teuchos::ParameterList plist;
  plist.set("Num Blocks",numBlocks_);
  plist.set("Recycled Blocks",recycleBlocks_);
  
  // Reset the status test.  
  outputTest_->reset();

  // Assume convergence is achieved, then let any failed convergence set this to false.
  bool isConverged = true;  

  // Compute AU = A*U, UTAU = U'*AU, AUTAU = (AU)'*(AU)
  if (existU_) {
    index.resize(recycleBlocks_);
    for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
    Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_,  index );
    index.resize(recycleBlocks_);
    for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
    Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
    // Initialize AU
    problem_->applyOp( *Utmp, *AUtmp );
    // Initialize UTAU
    Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
    MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
    // Initialize AUTAU  ( AUTAU = AU'*(M\AU) )
    Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
    if ( problem_->getLeftPrec() != Teuchos::null ) {
      index.resize(recycleBlocks_);
      for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
      index.resize(recycleBlocks_);
      for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
      Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage
      problem_->applyLeftPrec( *AUtmp, *LeftPCAU );
      MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
    } else {
      MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
    }
  }

  // RCG solver

  Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
  rcg_iter = Teuchos::rcp( new RCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );

  // Enter solve() iterations
  {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
    Teuchos::TimeMonitor slvtimer(*timerSolve_);
#endif

    while ( numRHS2Solve > 0 ) {

      // Debugging output to tell use if recycle space exists and will be used
      if (printer_->isVerbosity( Debug ) ) {
        if (existU_) printer_->print( Debug, "Using recycle space generated from previous call to solve()." );
        else printer_->print( Debug, "No recycle space exists." );
      }

      // Reset the number of iterations.
      rcg_iter->resetNumIters();

      // Set the current number of recycle blocks and subspace dimension with the RCG iteration.
      rcg_iter->setSize( recycleBlocks_, numBlocks_ );

      // Reset the number of calls that the status test output knows about.
      outputTest_->resetNumCalls();

      // indicate that updated recycle space has not yet been generated for this linear system
      existU1_ = false;

      // reset cycle count
      cycle = 0;

      // Get the current residual 
      problem_->computeCurrResVec( &*r_ );

      // If U exists, find best soln over this space first
      if (existU_) {
        // Solve linear system UTAU * y = (U'*r)
        Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
        index.resize(recycleBlocks_);
        for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
        Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_,  index );
        MVT::MvTransMv( one, *Utmp, *r_, Utr );
        Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
        Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
        LUUTAUtmp.assign(UTAUtmp);
        int info = 0;
        lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
        TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
                           "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");

        // Update solution (x = x + U*y)
        MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
 
        // Update residual ( r = r - AU*y )
        index.resize(recycleBlocks_);
        for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
        Teuchos::RCP<const MV> AUtmp  = MVT::CloneView( *AU_, index );
        MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
      }

      if ( problem_->getLeftPrec() != Teuchos::null ) {
        problem_->applyLeftPrec( *r_, *z_ );
      } else {
        z_ = r_;
      }
   
      // rTz_old = r'*z
      MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );

      if ( existU_ ) {
        // mu = UTAU\(AU'*z);
        Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
        index.resize(recycleBlocks_);
        for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
        Teuchos::RCP<const MV> AUtmp  = MVT::CloneView( *AU_, index );
        MVT::MvTransMv( one, *AUtmp, *z_, mu );
        Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
        char TRANS = 'N';
        int info;
        lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
        TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
                           "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
        // p  = z - U*mu;
        index.resize( 1 );
        index[0] = 0;
        Teuchos::RCP<MV> Ptmp  = MVT::CloneViewNonConst( *P_, index );
        MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
        MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
      } else {
        // p = z;
        index.resize( 1 );
        index[0] = 0;
        Teuchos::RCP<MV> Ptmp  = MVT::CloneViewNonConst( *P_, index );
        MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
      }

      // Set the new state and initialize the solver.
      RCGIterState<ScalarType,MV> newstate;

      // Create RCP views here 
      index.resize( numBlocks_+1 );
      for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
      newstate.P  = MVT::CloneViewNonConst( *P_,  index );
      index.resize( recycleBlocks_ );
      for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
      newstate.U  = MVT::CloneViewNonConst( *U_,  index );
      index.resize( recycleBlocks_ );
      for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
      newstate.AU  = MVT::CloneViewNonConst( *AU_,  index );
      newstate.Alpha = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) ); 
      newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) ); 
      newstate.D = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) ); 
      newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) ); 
      newstate.LUUTAU = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) ); 
      // assign the rest of the values in the struct
      newstate.curDim = 1; // We have initialized the first search vector
      newstate.Ap = Ap_;
      newstate.r = r_;
      newstate.z = z_;
      newstate.existU = existU_;
      newstate.ipiv = ipiv_;
      newstate.rTz_old = rTz_old_;

      rcg_iter->initialize(newstate);

      while(1) {

  // tell rcg_iter to iterate
  try {
    rcg_iter->iterate();

    //
    // check convergence first
    //
    if ( convTest_->getStatus() == Passed ) {
      // We have convergence
            break; // break from while(1){rcg_iter->iterate()}
    }
    //
    // check for maximum iterations
    //
    else if ( maxIterTest_->getStatus() == Passed ) {
      // we don't have convergence
      isConverged = false;
            break; // break from while(1){rcg_iter->iterate()}
    }
    //
    // check if cycle complete; update for next cycle
    //
    else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
            // index into P_ of last search vector generated this cycle
            int lastp = -1;
            // index into Beta_ of last entry generated this cycle
            int lastBeta = -1;
            if (recycleBlocks_ > 0) {
              if (!existU_) {
                if (cycle == 0) { // No U, no U1

                   Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
                   Ftmp.putScalar(zero);
                   Gtmp.putScalar(zero);
                   for (int ii=0;ii<numBlocks_;ii++) {
                     Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
                     if (ii > 0) {
                       Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
                       Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
                     }
                     Ftmp(ii,ii) = Dtmp(ii,0); 
                   }

                   // compute harmonic Ritz vectors
                   Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
                   getHarmonicVecs(Ftmp,Gtmp,Ytmp);

                   // U1 = [P(:,1:end-1)*Y];
                   index.resize( numBlocks_ );
                   for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_,  index );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
                   MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );

                   // Precompute some variables for next cycle

                   // AU1TAU1     = Y'*G*Y;
                   Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
                   GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
                   AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);

                   // AU1TU1      = Y'*F*Y;
                   Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
                   FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
                   AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);

                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
                   // Must reinitialize AU1TAP; can become dense later
                   AU1TAPtmp.putScalar(zero);
                   // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
                   ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
                   for (int ii=0; ii<recycleBlocks_; ++ii) {
                      AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
                   }
 
                   // indicate that updated recycle space now defined
                   existU1_ = true;

                   // Indicate the size of the P, Beta structures generated this cycle
                   lastp = numBlocks_;
                   lastBeta = numBlocks_-1;

                } // if (cycle == 0)
                else { // No U, but U1 guaranteed to exist now

                   // Finish computation of subblocks
                   // AU1TAP = AU1TAP * D(1);
                   Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
                   AU1TAPtmp.scale(Dtmp(0,0));

                   Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
                   APTAPtmp.putScalar(zero);
                   for (int ii=0; ii<numBlocks_; ii++) {
                     APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
                     if (ii > 0) {
                       APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
                       APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
                     }
                   }

                   // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
                   Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
                   F11.assign(AU1TU1tmp);
                   F12.putScalar(zero);
                   F21.putScalar(zero);
                   F22.putScalar(zero);
                   for(int ii=0;ii<numBlocks_;ii++) {
                     F22(ii,ii) = Dtmp(ii,0);
                   }

                   // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
                   Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
                   G11.assign(AU1TAU1tmp);
                   G12.assign(AU1TAPtmp);
                   // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
                   for (int ii=0;ii<recycleBlocks_;++ii)
                     for (int jj=0;jj<numBlocks_;++jj)
                       G21(jj,ii) = G12(ii,jj);
                   G22.assign(APTAPtmp);

                   // compute harmonic Ritz vectors
                   Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
                   getHarmonicVecs(Ftmp,Gtmp,Ytmp);

                   // U1 = [U1 P(:,2:end-1)]*Y;
                   index.resize( numBlocks_ );
                   for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
                   Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_,  index );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_,  index );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> U1Y1tmp  = MVT::CloneViewNonConst( *U1Y1_,  index );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
                   MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
                   MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
                   MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);

                   // Precompute some variables for next cycle
 
                   // AU1TAU1     = Y'*G*Y;
                   Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
                   GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
                   AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);

                   // AU1TAP      = zeros(k,m);
                   // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
                   AU1TAPtmp.putScalar(zero);
                   ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
                   for (int ii=0; ii<recycleBlocks_; ++ii) {
                      AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
                   }

                   // AU1TU1      = Y'*F*Y;
                   Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
                   FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
                   //Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
                   AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
 
                   // Indicate the size of the P, Beta structures generated this cycle
                   lastp = numBlocks_+1;
                   lastBeta = numBlocks_;

                } // if (cycle != 1)
              } // if (!existU_)
              else { // U exists
                if (cycle == 0) { // No U1, but U exists
                   Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
                   APTAPtmp.putScalar(zero);
                   for (int ii=0; ii<numBlocks_; ii++) {
                     APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
                     if (ii > 0) {
                       APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
                       APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
                     }
                   }

                   Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
                   L2tmp.putScalar(zero);
                   for(int ii=0;ii<numBlocks_;ii++) {
                     L2tmp(ii,ii)   = 1./Alphatmp(ii,0);
                     L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
                   }

                   // AUTAP = UTAU*Delta*L2;
                   Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
                   DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
                   AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);

                   // F = [UTAU zeros(k,m); zeros(m,k) diag(D)];
                   Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
                   F11.assign(UTAUtmp);
                   F12.putScalar(zero);
                   F21.putScalar(zero);
                   F22.putScalar(zero);
                   for(int ii=0;ii<numBlocks_;ii++) {
                     F22(ii,ii) = Dtmp(ii,0);
                   }
 
                   // G = [AUTAU AUTAP; AUTAP' APTAP];
                   Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
                   G11.assign(AUTAUtmp);
                   G12.assign(AUTAPtmp);
                   // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
                   for (int ii=0;ii<recycleBlocks_;++ii)
                     for (int jj=0;jj<numBlocks_;++jj)
                       G21(jj,ii) = G12(ii,jj);
                   G22.assign(APTAPtmp);

                   // compute harmonic Ritz vectors
                   Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
                   getHarmonicVecs(Ftmp,Gtmp,Ytmp);

                   // U1 = [U P(:,1:end-1)]*Y;
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
                   Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_,  index );
                   index.resize( numBlocks_ );
                   for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_,  index );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_,  index );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> UY1tmp  = MVT::CloneViewNonConst( *U1Y1_,  index );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
                   MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
                   MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
                   MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);

                   // Precompute some variables for next cycle
 
                   // AU1TAU1     = Y'*G*Y;
                   Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
                   GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
                   AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
 
                   // AU1TU1      = Y'*F*Y;
                   Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
                   FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
                   AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);

                   // AU1TU   = UTAU;
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
                   AU1TUtmp.assign(UTAUtmp);

                   // dold    = D(end);
                   dold = Dtmp(numBlocks_-1,0);

                   // indicate that updated recycle space now defined
                   existU1_ = true; 

                   // Indicate the size of the P, Beta structures generated this cycle
                   lastp = numBlocks_;
                   lastBeta = numBlocks_-1;
                }
                else { // Have U and U1
                   Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
                   for (int ii=0; ii<numBlocks_; ii++) {
                     APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
                     if (ii > 0) {
                       APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
                       APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
                     }
                   }

                   Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
                   for(int ii=0;ii<numBlocks_;ii++) {
                     L2tmp(ii,ii)   = 1./Alphatmp(ii,0);
                     L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
                   }

                   // M(end,1) = dold*(-Beta(1)/Alpha(1));
                   // AU1TAP = Y'*[AU1TU*Delta*L2; M];
                   Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
                   DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
                   AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
                   Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
                   AU1TAPtmp.putScalar(zero);
                   AU1TAPtmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
                   Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
                   for(int ii=0;ii<recycleBlocks_;ii++) {
                     AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
                   }

                   // AU1TU = Y1'*AU1TU
                   Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
                   Y1TAU1TU.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
                   AU1TUtmp.assign(Y1TAU1TU);

                   // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
                   Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
                   F11.assign(AU1TU1tmp);
                   for(int ii=0;ii<numBlocks_;ii++) {
                     F22(ii,ii) = Dtmp(ii,0);
                   }
 
                   // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
                   Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
                   Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
                   Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
                   G11.assign(AU1TAU1tmp);
                   G12.assign(AU1TAPtmp);
                   // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
                   for (int ii=0;ii<recycleBlocks_;++ii)
                     for (int jj=0;jj<numBlocks_;++jj)
                       G21(jj,ii) = G12(ii,jj);
                   G22.assign(APTAPtmp);

                   // compute harmonic Ritz vectors
                   Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
                   getHarmonicVecs(Ftmp,Gtmp,Ytmp);

                   // U1 = [U1 P(:,2:end-1)]*Y;
                   index.resize( numBlocks_ );
                   for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
                   Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_,  index );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_,  index );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> U1tmp  = MVT::CloneViewNonConst( *U1_,  index );
                   index.resize( recycleBlocks_ );
                   for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
                   Teuchos::RCP<MV> U1Y1tmp  = MVT::CloneViewNonConst( *U1Y1_,  index );
                   MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
                   MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
                   MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
 
                   // Precompute some variables for next cycle

                   // AU1TAU1     = Y'*G*Y;
                   Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
                   GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
                   AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
 
                   // AU1TU1      = Y'*F*Y;
                   Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
                   FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
                   AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);

                   // dold    = D(end);
                   dold = Dtmp(numBlocks_-1,0);

                   // Indicate the size of the P, Beta structures generated this cycle
                   lastp = numBlocks_+1;
                   lastBeta = numBlocks_;

                }
              }
            } // if (recycleBlocks_ > 0)

            // Cleanup after end of cycle

            // P = P(:,end-1:end);
            index.resize( 2 ); 
            index[0] = lastp-1; index[1] = lastp;
            Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_,  index );
            index[0] = 0;       index[1] = 1;
            MVT::SetBlock(*Ptmp2,index,*P_);

            // Beta = Beta(end);
            (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
 
            // Delta = Delta(:,end);
            if (existU_) { // Delta only initialized if U exists
              Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
              Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
              mu1.assign(mu2);
            }

            // Now reinitialize state variables for next cycle
            newstate.P = Teuchos::null;
            index.resize( numBlocks_+1 );
            for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
            newstate.P  = MVT::CloneViewNonConst( *P_,  index );

            newstate.Beta = Teuchos::null;
            newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );

            newstate.Delta = Teuchos::null;
            newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );

            newstate.curDim = 1; // We have initialized the first search vector

            // Pass to iteration object
            rcg_iter->initialize(newstate);

            // increment cycle count
            cycle = cycle + 1;

    }
    //
    // we returned from iterate(), but none of our status tests Passed.
    // something is wrong, and it is probably our fault.
    //
    else {
      TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
             "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
    }
  }
  catch (const std::exception &e) {
    printer_->stream(Errors) << "Error! Caught std::exception in RCGIter::iterate() at iteration " 
           << rcg_iter->getNumIters() << std::endl 
           << e.what() << std::endl;
    throw;
  }
      }

      // Inform the linear problem that we are finished with this block linear system.
      problem_->setCurrLS();
      
      // Update indices for the linear systems to be solved.
      numRHS2Solve -= 1;
      if ( numRHS2Solve > 0 ) {
        currIdx[0]++;
  // Set the next indices.
  problem_->setLSIndex( currIdx );
      } 
      else {
        currIdx.resize( numRHS2Solve );
      }

      // Update the recycle space for the next linear system
      if (existU1_) { // be sure updated recycle space was created
        // U = U1
        index.resize(recycleBlocks_);
        for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
        MVT::SetBlock(*U1_,index,*U_);
        // Set flag indicating recycle space is now defined
        existU_ = true;
        if (numRHS2Solve > 0) { // also update AU, UTAU, and AUTAU
          // Free pointers in newstate
          newstate.P = Teuchos::null;
          newstate.Ap = Teuchos::null;
          newstate.r = Teuchos::null;
          newstate.z = Teuchos::null;
          newstate.U = Teuchos::null;
          newstate.AU = Teuchos::null;
          newstate.Alpha = Teuchos::null;
          newstate.Beta = Teuchos::null;
          newstate.D = Teuchos::null;
          newstate.Delta = Teuchos::null;
          newstate.LUUTAU = Teuchos::null;
          newstate.ipiv = Teuchos::null;
          newstate.rTz_old = Teuchos::null;

          // Reinitialize AU, UTAU, AUTAU
          index.resize(recycleBlocks_);
          for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
          Teuchos::RCP<const MV> Utmp  = MVT::CloneView( *U_,  index );
          index.resize(recycleBlocks_);
          for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
          Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
          // Initialize AU
          problem_->applyOp( *Utmp, *AUtmp );
          // Initialize UTAU
          Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
          MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
          // Initialize AUTAU  ( AUTAU = AU'*(M\AU) )
          Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
          if ( problem_->getLeftPrec() != Teuchos::null ) {
            index.resize(recycleBlocks_);
            for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
            index.resize(recycleBlocks_);
            for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
            Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage
            problem_->applyLeftPrec( *AUtmp, *LeftPCAU );
            MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
          } else {
            MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
          }
        } // if (numRHS2Solve > 0)

      } // if (existU1) 
    }// while ( numRHS2Solve > 0 )
    
  }

  // print final summary
  sTest_->print( printer_->stream(FinalSummary) );
 
  // print timing information
#ifdef BELOS_TEUCHOS_TIME_MONITOR
  // Calling summarize() can be expensive, so don't call unless the
  // user wants to print out timing details.  summarize() will do all
  // the work even if it's passed a "black hole" output stream.
  if (verbosity_ & TimingDetails)
    Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
#endif

  // get iteration information for this solve
  numIters_ = maxIterTest_->getNumIters();

  // Save the convergence test value ("achieved tolerance") for this solve.
  {
    using Teuchos::rcp_dynamic_cast;
    typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
    // testValues is nonnull and not persistent.
    const std::vector<MagnitudeType>* pTestValues = 
      rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
    
    TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
      "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
      "method returned NULL.  Please report this bug to the Belos developers.");
    
    TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
      "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
      "method returned a vector of length zero.  Please report this bug to the "
      "Belos developers.");

    // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
    // achieved tolerances for all vectors in the current solve(), or
    // just for the vectors from the last deflation?
    achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
  }
  
  if (!isConverged) {
    return Unconverged; // return from RCGSolMgr::solve() 
  }
  return Converged; // return from RCGSolMgr::solve() 
}

//  Compute the harmonic eigenpairs of the projected, dense system.
template<class ScalarType, class MV, class OP>
void RCGSolMgr<ScalarType,MV,OP>::getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType>& F,
                                                  const Teuchos::SerialDenseMatrix<int,ScalarType>& G,
                                                        Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
  // order of F,G 
  int n = F.numCols(); 

  // The LAPACK interface
  Teuchos::LAPACK<int,ScalarType> lapack;
 
  // Magnitude of harmonic Ritz values
  std::vector<MagnitudeType> w(n);
 
  // Sorted order of harmonic Ritz values
  std::vector<int> iperm(n);
 
  // Compute k smallest harmonic Ritz pairs
  // SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO )
  int itype = 1; // solve A*x = (lambda)*B*x
  char jobz='V'; // compute eigenvalues and eigenvectors
  char uplo='U'; // since F,G symmetric, reference only their upper triangular data
  std::vector<ScalarType> work(1);
  int lwork = -1;
  int info = 0;
  // since SYGV destroys workspace, create copies of F,G
  Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_, F_->numRows(), F_->numCols() );
  Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_, G_->numRows(), G_->numCols() );

  // query for optimal workspace size
  lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info); 
  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
                     "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
  lwork = (int)work[0];
  work.resize(lwork);
  lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info); 
  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
                     "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");


  // Construct magnitude of each harmonic Ritz value
  this->sort(w,n,iperm);
 
  // Select recycledBlocks_ smallest eigenvectors
  for( int i=0; i<recycleBlocks_; i++ ) {
    for( int j=0; j<n; j++ ) {
      Y(j,i) = G2(j,iperm[i]);
    }
  }
 
}

// This method sorts list of n floating-point numbers and return permutation vector
template<class ScalarType, class MV, class OP>
void RCGSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm)
{
  int l, r, j, i, flag;
  int    RR2;
  double dRR, dK;
 
  // Initialize the permutation vector.
  for(j=0;j<n;j++)
    iperm[j] = j;
 
  if (n <= 1) return;
 
  l    = n / 2 + 1;
  r    = n - 1;
  l    = l - 1;
  dRR  = dlist[l - 1];
  dK   = dlist[l - 1];
 
  RR2 = iperm[l - 1];
  while (r != 0) {
    j = l;
    flag = 1;
 
    while (flag == 1) {
      i = j;
      j = j + j;
 
      if (j > r + 1)
        flag = 0;
      else {
        if (j < r + 1)
          if (dlist[j] > dlist[j - 1]) j = j + 1;
 
        if (dlist[j - 1] > dK) {
          dlist[i - 1] = dlist[j - 1];
          iperm[i - 1] = iperm[j - 1];
        }
        else {
          flag = 0;
        }
      }
    }
    dlist[i - 1] = dRR;
    iperm[i - 1] = RR2;
    if (l == 1) {
      dRR  = dlist [r];
      RR2 = iperm[r];
      dK = dlist[r];
      dlist[r] = dlist[0];
      iperm[r] = iperm[0];
      r = r - 1;
    }
    else {
      l   = l - 1;
      dRR  = dlist[l - 1];
      RR2  = iperm[l - 1];
      dK   = dlist[l - 1];
    }
  }
  dlist[0] = dRR;
  iperm[0] = RR2;
}

//  This method requires the solver manager to return a std::string that describes itself.
template<class ScalarType, class MV, class OP>
std::string RCGSolMgr<ScalarType,MV,OP>::description() const
{
  std::ostringstream oss;
  oss << "Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
  return oss.str();
}
  
} // end Belos namespace

#endif /* BELOS_RCG_SOLMGR_HPP */