AnasaziLOBPCG.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//                 Anasazi: Block Eigensolvers Package
//
// Copyright 2004 NTESS and the Anasazi contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER


/*
    LOBPCG contains local storage of up to 10*blockSize_ vectors, representing 10 entities
      X,H,P,R
      KX,KH,KP  (product of K and the above)
      MX,MH,MP  (product of M and the above, not allocated if we don't have an M matrix)
    If full orthogonalization is enabled, one extra multivector of blockSize_ vectors is required to 
    compute the local update of X and P.
    
    A solver is bound to an eigenproblem at declaration.
    Other solver parameters (e.g., block size, auxiliary vectors) can be changed dynamically.
    
    The orthogonalization manager is used to project away from the auxiliary vectors.
    If full orthogonalization is enabled, the orthogonalization manager is also used to construct an M orthonormal basis.
    The orthogonalization manager is subclass of MatOrthoManager, which LOBPCG assumes to be defined by the M inner product.
    LOBPCG will not work correctly if the orthomanager uses a different inner product.
 */


#ifndef ANASAZI_LOBPCG_HPP
#define ANASAZI_LOBPCG_HPP

#include "AnasaziTypes.hpp"

#include "AnasaziEigensolver.hpp"
#include "AnasaziMultiVecTraits.hpp"
#include "AnasaziOperatorTraits.hpp"
#include "Teuchos_ScalarTraits.hpp"

#include "AnasaziMatOrthoManager.hpp"
#include "AnasaziSolverUtils.hpp"

#include "Teuchos_LAPACK.hpp"
#include "Teuchos_BLAS.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_TimeMonitor.hpp"

namespace Anasazi {



  template <class ScalarType, class MultiVector>
  struct LOBPCGState {
    Teuchos::RCP<const MultiVector> V; 
    Teuchos::RCP<const MultiVector> KV; 
    Teuchos::RCP<const MultiVector> MV;

    Teuchos::RCP<const MultiVector> X; 
    Teuchos::RCP<const MultiVector> KX; 
    Teuchos::RCP<const MultiVector> MX;

    Teuchos::RCP<const MultiVector> P; 
    Teuchos::RCP<const MultiVector> KP; 
    Teuchos::RCP<const MultiVector> MP;

    Teuchos::RCP<const MultiVector> H; 
    Teuchos::RCP<const MultiVector> KH; 
    Teuchos::RCP<const MultiVector> MH;

    Teuchos::RCP<const MultiVector> R;

    Teuchos::RCP<const std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > T;

    LOBPCGState() : 
      V(Teuchos::null),KV(Teuchos::null),MV(Teuchos::null),
      X(Teuchos::null),KX(Teuchos::null),MX(Teuchos::null),
      P(Teuchos::null),KP(Teuchos::null),MP(Teuchos::null),
      H(Teuchos::null),KH(Teuchos::null),MH(Teuchos::null),
      R(Teuchos::null),T(Teuchos::null) {};
  };




  class LOBPCGRitzFailure : public AnasaziError {public:
    LOBPCGRitzFailure(const std::string& what_arg) : AnasaziError(what_arg)
    {}};

  class LOBPCGInitFailure : public AnasaziError {public:
    LOBPCGInitFailure(const std::string& what_arg) : AnasaziError(what_arg)
    {}};

  class LOBPCGOrthoFailure : public AnasaziError {public:
    LOBPCGOrthoFailure(const std::string& what_arg) : AnasaziError(what_arg)
    {}};



  template <class ScalarType, class MV, class OP>
  class LOBPCG : public Eigensolver<ScalarType,MV,OP> { 
  public:
    

    
    LOBPCG( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 
            const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
            const Teuchos::RCP<OutputManager<ScalarType> > &printer,
            const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
            const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
            Teuchos::ParameterList &params 
          );
    
    virtual ~LOBPCG() {};




    void iterate();

    void initialize(LOBPCGState<ScalarType,MV>& newstate);

    void initialize();

    bool isInitialized() const;

    LOBPCGState<ScalarType,MV> getState() const;




    int getNumIters() const;

    void resetNumIters();

    Teuchos::RCP<const MV> getRitzVectors();

    std::vector<Value<ScalarType> > getRitzValues();

    std::vector<int> getRitzIndex();


    std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getResNorms();


    std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRes2Norms();


    std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> getRitzRes2Norms();


    int getCurSubspaceDim() const;

    int getMaxSubspaceDim() const;




    void setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test);

    Teuchos::RCP<StatusTest<ScalarType,MV,OP> > getStatusTest() const;

    const Eigenproblem<ScalarType,MV,OP>& getProblem() const;


    void setBlockSize(int blockSize);


    int getBlockSize() const;


    void setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs);

    Teuchos::Array<Teuchos::RCP<const MV> > getAuxVecs() const;




    void setFullOrtho(bool fullOrtho);

    bool getFullOrtho() const;
    
    bool hasP();

    


    void currentStatus(std::ostream &os);


  private:
    //
    //
    //
    void setupViews();
    //
    // Convenience typedefs
    //
    typedef SolverUtils<ScalarType,MV,OP> Utils;
    typedef MultiVecTraits<ScalarType,MV> MVT;
    typedef OperatorTraits<ScalarType,MV,OP> OPT;
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef typename SCT::magnitudeType MagnitudeType;
    const MagnitudeType ONE;  
    const MagnitudeType ZERO; 
    const MagnitudeType NANVAL;
    //
    // Internal structs
    //
    struct CheckList {
      bool checkX, checkMX, checkKX;
      bool checkH, checkMH;
      bool checkP, checkMP, checkKP;
      bool checkR, checkQ;
      CheckList() : checkX(false),checkMX(false),checkKX(false),
                    checkH(false),checkMH(false),
                    checkP(false),checkMP(false),checkKP(false),
                    checkR(false),checkQ(false) {};
    };
    //
    // Internal methods
    //
    std::string accuracyCheck(const CheckList &chk, const std::string &where) const;
    //
    // Classes inputed through constructor that define the eigenproblem to be solved.
    //
    const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> >     problem_;
    const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > sm_;
    const Teuchos::RCP<OutputManager<ScalarType> >          om_;
    Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       tester_;
    const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >  orthman_;
    //
    // Information obtained from the eigenproblem
    //
    Teuchos::RCP<const OP> Op_;
    Teuchos::RCP<const OP> MOp_;
    Teuchos::RCP<const OP> Prec_;
    bool hasM_;
    //
    // Internal timers
    //
    Teuchos::RCP<Teuchos::Time> timerOp_, timerMOp_, timerPrec_,
                                timerSort_, 
                                timerLocalProj_, timerDS_,
                                timerLocalUpdate_, timerCompRes_,
                                timerOrtho_, timerInit_;
    //
    // Counters
    //
    // Number of operator applications
    int count_ApplyOp_, count_ApplyM_, count_ApplyPrec_;

    //
    // Algorithmic parameters.
    //
    // blockSize_ is the solver block size
    int blockSize_;
    //
    // fullOrtho_ dictates whether the orthogonalization procedures specified by Hetmaniuk and Lehoucq should
    // be activated (see citations at the top of this file)
    bool fullOrtho_;

    //
    // Current solver state
    //
    // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
    // is capable of running; _initialize is controlled  by the initialize() member method
    // For the implications of the state of initialized_, please see documentation for initialize()
    bool initialized_;
    //
    // nevLocal_ reflects how much of the current basis is valid (0 <= nevLocal_ <= 3*blockSize_)
    // this tells us how many of the values in theta_ are valid Ritz values
    int nevLocal_;
    //
    // hasP_ tells us whether there is valid data in P (and KP,MP)
    bool hasP_;
    //
    // State Multivecs
    // V_, KV_ MV_  and  R_ are primary pointers to allocated multivectors
    Teuchos::RCP<MV> V_, KV_, MV_, R_;
    // the rest are multivector views into V_, KV_ and MV_
    Teuchos::RCP<MV> X_, KX_, MX_,
                     H_, KH_, MH_,
                     P_, KP_, MP_;

    //
    // if fullOrtho_ == true, then we must produce the following on every iteration:
    // [newX newP] = [X H P] [CX;CP]
    // the structure of [CX;CP] when using full orthogonalization does not allow us to 
    // do this in situ, and R_ does not have enough storage for newX and newP. therefore, 
    // we must allocate additional storage for this.
    // otherwise, when not using full orthogonalization, the structure
    // [newX newP] = [X H P] [CX1  0 ]
    //                       [CX2 CP2]  allows us to work using only R as work space
    //                       [CX3 CP3] 
    Teuchos::RCP<MV> tmpmvec_;        
    // 
    // auxiliary vectors
    Teuchos::Array<Teuchos::RCP<const MV> > auxVecs_;
    int numAuxVecs_;
    //
    // Number of iterations that have been performed.
    int iter_;
    // 
    // Current eigenvalues, residual norms
    std::vector<MagnitudeType> theta_, Rnorms_, R2norms_;
    // 
    // are the residual norms current with the residual?
    bool Rnorms_current_, R2norms_current_;

  };




  // Constructor
  template <class ScalarType, class MV, class OP>
  LOBPCG<ScalarType,MV,OP>::LOBPCG(
        const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, 
        const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
        const Teuchos::RCP<OutputManager<ScalarType> > &printer,
        const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
        const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
        Teuchos::ParameterList &params
        ) :
    ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
    ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
    NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan()),
    // problem, tools
    problem_(problem), 
    sm_(sorter),
    om_(printer),
    tester_(tester),
    orthman_(ortho),
    // timers, counters
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
    timerOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Operation Op*x")),
    timerMOp_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Operation M*x")),
    timerPrec_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Operation Prec*x")),
    timerSort_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Sorting eigenvalues")),
    timerLocalProj_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Local projection")),
    timerDS_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Direct solve")),
    timerLocalUpdate_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Local update")),
    timerCompRes_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Computing residuals")),
    timerOrtho_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Orthogonalization")),
    timerInit_(Teuchos::TimeMonitor::getNewTimer("Anasazi: LOBPCG::Initialization")),
#endif
    count_ApplyOp_(0),
    count_ApplyM_(0),
    count_ApplyPrec_(0),
    // internal data
    blockSize_(0),
    fullOrtho_(params.get("Full Ortho", true)),
    initialized_(false),
    nevLocal_(0),
    hasP_(false),
    auxVecs_( Teuchos::Array<Teuchos::RCP<const MV> >(0) ), 
    numAuxVecs_(0),
    iter_(0),
    Rnorms_current_(false),
    R2norms_current_(false)
  {     
    TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,std::invalid_argument,
                       "Anasazi::LOBPCG::constructor: user passed null problem pointer.");
    TEUCHOS_TEST_FOR_EXCEPTION(sm_ == Teuchos::null,std::invalid_argument,
                       "Anasazi::LOBPCG::constructor: user passed null sort manager pointer.");
    TEUCHOS_TEST_FOR_EXCEPTION(om_ == Teuchos::null,std::invalid_argument,
                       "Anasazi::LOBPCG::constructor: user passed null output manager pointer.");
    TEUCHOS_TEST_FOR_EXCEPTION(tester_ == Teuchos::null,std::invalid_argument,
                       "Anasazi::LOBPCG::constructor: user passed null status test pointer.");
    TEUCHOS_TEST_FOR_EXCEPTION(orthman_ == Teuchos::null,std::invalid_argument,
                       "Anasazi::LOBPCG::constructor: user passed null orthogonalization manager pointer.");
    TEUCHOS_TEST_FOR_EXCEPTION(problem_->isProblemSet() == false, std::invalid_argument,
                       "Anasazi::LOBPCG::constructor: problem is not set.");
    TEUCHOS_TEST_FOR_EXCEPTION(problem_->isHermitian() == false, std::invalid_argument,
                       "Anasazi::LOBPCG::constructor: problem is not Hermitian; LOBPCG requires Hermitian problem.");

    // get the problem operators
    Op_   = problem_->getOperator();
    TEUCHOS_TEST_FOR_EXCEPTION(Op_ == Teuchos::null, std::invalid_argument,
                       "Anasazi::LOBPCG::constructor: problem provides no operator.");
    MOp_  = problem_->getM();
    Prec_ = problem_->getPrec();
    hasM_ = (MOp_ != Teuchos::null);

    // set the block size and allocate data
    int bs = params.get("Block Size", problem_->getNEV());
    setBlockSize(bs);
  }


  // Set the block size and make necessary adjustments.
  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::setBlockSize (int newBS) 
  {
    // time spent here counts towards timerInit_
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
    Teuchos::TimeMonitor lcltimer( *timerInit_ );
#endif

    // This routine only allocates space; it doesn't not perform any computation
    // if size is decreased, take the first newBS vectors of all and leave state as is
    // otherwise, grow/allocate space and set solver to unitialized

    Teuchos::RCP<const MV> tmp;
    // grab some Multivector to Clone
    // in practice, getInitVec() should always provide this, but it is possible to use a 
    // Eigenproblem with nothing in getInitVec() by manually initializing with initialize(); 
    // in case of that strange scenario, we will try to Clone from R_ because it is smaller 
    // than V_, and we don't want to keep V_ around longer than necessary
    if (blockSize_ > 0) {
      tmp = R_;
    }
    else {
      tmp = problem_->getInitVec();
      TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::logic_error,
                         "Anasazi::LOBPCG::setBlockSize(): eigenproblem did not specify initial vectors to clone from.");
    }
    
    TEUCHOS_TEST_FOR_EXCEPTION(newBS <= 0 || static_cast<ptrdiff_t>(newBS) > MVT::GetGlobalLength(*tmp), 
                         std::invalid_argument, "Anasazi::LOBPCG::setBlockSize(): block size must be strictly positive.");
    if (newBS == blockSize_) {
      // do nothing
      return;
    }
    else if (newBS < blockSize_ && initialized_) {
      //
      // shrink vectors

      // release views so we can modify the bases
      X_ = Teuchos::null;
      KX_ = Teuchos::null;
      MX_ = Teuchos::null;
      H_ = Teuchos::null;
      KH_ = Teuchos::null;
      MH_ = Teuchos::null;
      P_ = Teuchos::null;
      KP_ = Teuchos::null;
      MP_ = Teuchos::null;

      // make new indices vectors
      std::vector<int> newind(newBS), oldind(newBS);
      for (int i=0; i<newBS; i++) {
        newind[i] = i;
        oldind[i] = i;
      }

      Teuchos::RCP<MV> newV, newMV, newKV, newR;
      Teuchos::RCP<const MV> src;
      // allocate R and newV
      newR = MVT::Clone(*tmp,newBS);
      newV = MVT::Clone(*tmp,newBS*3);
      newKV = MVT::Clone(*tmp,newBS*3);
      if (hasM_) {
        newMV = MVT::Clone(*tmp,newBS*3);
      }

      //
      // if we are initialized, we want to pull the data from V_ into newV:
      //           bs  |  bs  | bs 
      // newV  = [newX | **** |newP ]
      // newKV = [newKX| **** |newKP]
      // newMV = [newMX| **** |newMP]
      // where 
      //          oldbs   |  oldbs  |   oldbs   
      //  V_ = [newX  *** | ******* | newP  ***]
      // KV_ = [newKX *** | ******* | newKP ***]
      // MV_ = [newMX *** | ******* | newMP ***]
      //
      // we don't care to copy the data corresponding to H
      // we will not copy the M data if !hasM_, because it doesn't exist
      //

      // these are shrink operations which preserve their data
      theta_.resize(3*newBS);
      Rnorms_.resize(newBS);
      R2norms_.resize(newBS);

      // copy residual vectors: oldind,newind currently contains [0,...,newBS-1]
      src = MVT::CloneView(*R_,newind);
      MVT::SetBlock(*src,newind,*newR);
      // free old memory and point to new memory
      R_ = newR;

      // copy in order: newX newKX newMX, then newP newKP newMP
      // for  X: [0,bs-1] <-- [0,bs-1] 
      src = MVT::CloneView(*V_,oldind);
      MVT::SetBlock(*src,newind,*newV);
      src = MVT::CloneView(*KV_,oldind);
      MVT::SetBlock(*src,newind,*newKV);
      if (hasM_) {
        src = MVT::CloneView(*MV_,oldind);
        MVT::SetBlock(*src,newind,*newMV);
      }
      // for  P: [2*bs, 3*bs-1] <-- [2*oldbs, 2*oldbs+bs-1] 
      for (int i=0; i<newBS; i++) {
        newind[i] += 2*newBS;
        oldind[i] += 2*blockSize_;
      }
      src = MVT::CloneView(*V_,oldind);
      MVT::SetBlock(*src,newind,*newV);
      src = MVT::CloneView(*KV_,oldind);
      MVT::SetBlock(*src,newind,*newKV);
      if (hasM_) {
        src = MVT::CloneView(*MV_,oldind);
        MVT::SetBlock(*src,newind,*newMV);
      }

      // release temp view
      src = Teuchos::null;

      // release old allocations and point at new memory
      V_ = newV;
      KV_ = newKV;
      if (hasM_) {
        MV_ = newMV;
      }
      else {
        MV_ = V_;
      }
    }
    else {  
      // newBS > blockSize_  or  not initialized
      // this is also the scenario for our initial call to setBlockSize(), in the constructor
      initialized_ = false;
      hasP_ = false;

      // release views
      X_ = Teuchos::null;
      KX_ = Teuchos::null;
      MX_ = Teuchos::null;
      H_ = Teuchos::null;
      KH_ = Teuchos::null;
      MH_ = Teuchos::null;
      P_ = Teuchos::null;
      KP_ = Teuchos::null;
      MP_ = Teuchos::null;

      // free allocated storage
      R_ = Teuchos::null;
      V_ = Teuchos::null;

      // allocate scalar vectors
      theta_.resize(3*newBS,NANVAL);
      Rnorms_.resize(newBS,NANVAL);
      R2norms_.resize(newBS,NANVAL);
      
      // clone multivectors off of tmp
      R_ = MVT::Clone(*tmp,newBS);
      V_ = MVT::Clone(*tmp,newBS*3);
      KV_ = MVT::Clone(*tmp,newBS*3);
      if (hasM_) {
        MV_ = MVT::Clone(*tmp,newBS*3);
      }
      else {
        MV_ = V_;
      }
    }

    // allocate tmp space
    tmpmvec_ = Teuchos::null;
    if (fullOrtho_) {
      tmpmvec_ = MVT::Clone(*tmp,newBS);
    }

    // set new block size
    blockSize_ = newBS;

    // setup new views
    setupViews();
  }


  // Setup views into V,KV,MV
  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::setupViews() 
  {
    std::vector<int> ind(blockSize_);

    for (int i=0; i<blockSize_; i++) {
      ind[i] = i;
    }
    X_  = MVT::CloneViewNonConst(*V_,ind);
    KX_ = MVT::CloneViewNonConst(*KV_,ind);
    if (hasM_) {
      MX_ = MVT::CloneViewNonConst(*MV_,ind);
    }
    else {
      MX_ = X_;
    }

    for (int i=0; i<blockSize_; i++) {
      ind[i] += blockSize_;
    }
    H_  = MVT::CloneViewNonConst(*V_,ind);
    KH_ = MVT::CloneViewNonConst(*KV_,ind);
    if (hasM_) {
      MH_ = MVT::CloneViewNonConst(*MV_,ind);
    }
    else {
      MH_ = H_;
    }

    for (int i=0; i<blockSize_; i++) {
      ind[i] += blockSize_;
    }
    P_  = MVT::CloneViewNonConst(*V_,ind);
    KP_ = MVT::CloneViewNonConst(*KV_,ind);
    if (hasM_) {
      MP_ = MVT::CloneViewNonConst(*MV_,ind);
    }
    else {
      MP_ = P_;
    }
  }


  // Set the auxiliary vectors
  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::setAuxVecs(const Teuchos::Array<Teuchos::RCP<const MV> > &auxvecs) {
    typedef typename Teuchos::Array<Teuchos::RCP<const MV> >::iterator tarcpmv;

    // set new auxiliary vectors
    auxVecs_ = auxvecs;
    
    numAuxVecs_ = 0;
    for (tarcpmv i=auxVecs_.begin(); i != auxVecs_.end(); i++) {
      numAuxVecs_ += MVT::GetNumberVecs(**i);
    }
    
    // If the solver has been initialized, X and P are not necessarily orthogonal to new auxiliary vectors
    if (numAuxVecs_ > 0 && initialized_) {
      initialized_ = false;
      hasP_ = false;
    }

    if (om_->isVerbosity( Debug ) ) {
      // Check almost everything here
      CheckList chk;
      chk.checkQ = true;
      om_->print( Debug, accuracyCheck(chk, ": in setAuxVecs()") );
    }
  }


  /* Initialize the state of the solver
   * 
   * POST-CONDITIONS:
   *
   * initialized_ == true
   * X is orthonormal, orthogonal to auxVecs_
   * KX = Op*X
   * MX = M*X if hasM_
   * theta_ contains Ritz values of X
   * R = KX - MX*diag(theta_)
   * if hasP() == true,
   *   P orthogonal to auxVecs_
   *   if fullOrtho_ == true,
   *     P orthonormal and orthogonal to X
   *   KP = Op*P
   *   MP = M*P
   */
  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::initialize(LOBPCGState<ScalarType,MV>& newstate)
  {
    // NOTE: memory has been allocated by setBlockSize(). Use SetBlock below; do not Clone
    // NOTE: Overall time spent in this routine is counted to timerInit_; portions will also be counted towards other primitives

#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
    Teuchos::TimeMonitor inittimer( *timerInit_ );
#endif

    std::vector<int> bsind(blockSize_);
    for (int i=0; i<blockSize_; i++) bsind[i] = i;

    // in LOBPCG, X (the subspace iterate) is primary
    // the order of dependence follows like so.
    // --init->                 X
    //    --op apply->          MX,KX
    //       --ritz analysis->  theta
    //          --optional->    P,MP,KP
    // 
    // if the user specifies all data for a level, we will accept it.
    // otherwise, we will generate the whole level, and all subsequent levels.
    //
    // the data members are ordered based on dependence, and the levels are
    // partitioned according to the amount of work required to produce the
    // items in a level.
    //
    // inconsitent multivectors widths and lengths will not be tolerated, and
    // will be treated with exceptions.

    // set up X, KX, MX: get them from "state" if user specified them

    //----------------------------------------
    // set up X, MX, KX
    //----------------------------------------
    if (newstate.X != Teuchos::null) {
      TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.X) != MVT::GetGlobalLength(*X_),
                          std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.X not correct." );
      // newstate.X must have blockSize_ vectors; any more will be ignored
      TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.X) < blockSize_,
                          std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.X must have at least block size vectors.");

      // put X data in X_
      MVT::SetBlock(*newstate.X,bsind,*X_);

      // put MX data in MX_
      if (hasM_) {
        if (newstate.MX != Teuchos::null) {
          TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MX) != MVT::GetGlobalLength(*MX_),
                              std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.MX not correct." );
          // newstate.MX must have blockSize_ vectors; any more will be ignored
          TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MX) < blockSize_,
                              std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.MX must have at least block size vectors.");
          MVT::SetBlock(*newstate.MX,bsind,*MX_);
        }
        else {
          // user didn't specify MX, compute it
          {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
            Teuchos::TimeMonitor lcltimer( *timerMOp_ );
#endif
            OPT::Apply(*MOp_,*X_,*MX_);
            count_ApplyM_ += blockSize_;
          }
          // we generated MX; we will generate R as well
          newstate.R = Teuchos::null;
        }
      }
  
      // put data in KX
      if (newstate.KX != Teuchos::null) {
        TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KX) != MVT::GetGlobalLength(*KX_),
                            std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.KX not correct." );
        // newstate.KX must have blockSize_ vectors; any more will be ignored
        TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KX) < blockSize_,
                            std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.KX must have at least block size vectors.");
        MVT::SetBlock(*newstate.KX,bsind,*KX_);
      }
      else {
        // user didn't specify KX, compute it
        {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
          Teuchos::TimeMonitor lcltimer( *timerOp_ );
#endif
          OPT::Apply(*Op_,*X_,*KX_);
          count_ApplyOp_ += blockSize_;
        }
        // we generated KX; we will generate R as well
        newstate.R = Teuchos::null;
      }
    }
    else {
      // user did not specify X
      // we will initialize X, compute KX and MX, and compute R
      //
      // clear state so we won't use any data from it below
      newstate.P  = Teuchos::null;
      newstate.KP = Teuchos::null;
      newstate.MP = Teuchos::null;
      newstate.R  = Teuchos::null;
      newstate.T  = Teuchos::null;

      // generate a basis and projectAndNormalize
      Teuchos::RCP<const MV> ivec = problem_->getInitVec();
      TEUCHOS_TEST_FOR_EXCEPTION(ivec == Teuchos::null,std::logic_error,
                         "Anasazi::LOBPCG::initialize(): Eigenproblem did not specify initial vectors to clone from.");

      int initSize = MVT::GetNumberVecs(*ivec);
      if (initSize > blockSize_) {
        // we need only the first blockSize_ vectors from ivec; get a view of them
        initSize = blockSize_;
        std::vector<int> ind(blockSize_);
        for (int i=0; i<blockSize_; i++) ind[i] = i;
        ivec = MVT::CloneView(*ivec,ind);
      }

      // assign ivec to first part of X_ 
      if (initSize > 0) {
        std::vector<int> ind(initSize);
        for (int i=0; i<initSize; i++) ind[i] = i;
        MVT::SetBlock(*ivec,ind,*X_);
      }
      // fill the rest of X_ with random
      if (blockSize_ > initSize) {
        std::vector<int> ind(blockSize_ - initSize);
        for (int i=0; i<blockSize_ - initSize; i++) ind[i] = initSize + i;
        Teuchos::RCP<MV> rX = MVT::CloneViewNonConst(*X_,ind);
        MVT::MvRandom(*rX);
        rX = Teuchos::null;
      }

      // put data in MX
      if (hasM_) {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerMOp_ );
#endif
        OPT::Apply(*MOp_,*X_,*MX_);
        count_ApplyM_ += blockSize_;
      }
  
      // remove auxVecs from X_ and normalize it
      if (numAuxVecs_ > 0) {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
#endif
        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummy;
        int rank = orthman_->projectAndNormalizeMat(*X_,auxVecs_,dummy,Teuchos::null,MX_);
        TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_, LOBPCGInitFailure,
                           "Anasazi::LOBPCG::initialize(): Couldn't generate initial basis of full rank.");
      }
      else {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
#endif
        int rank = orthman_->normalizeMat(*X_,Teuchos::null,MX_);
        TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_, LOBPCGInitFailure,
                           "Anasazi::LOBPCG::initialize(): Couldn't generate initial basis of full rank.");
      }

      // put data in KX
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerOp_ );
#endif
        OPT::Apply(*Op_,*X_,*KX_);
        count_ApplyOp_ += blockSize_;
      }
    } // end if (newstate.X != Teuchos::null)


    //----------------------------------------
    // set up Ritz values
    //----------------------------------------
    theta_.resize(3*blockSize_,NANVAL);
    if (newstate.T != Teuchos::null) {
      TEUCHOS_TEST_FOR_EXCEPTION( (signed int)(newstate.T->size()) < blockSize_,
                          std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.T must contain at least block size Ritz values.");
      for (int i=0; i<blockSize_; i++) {
        theta_[i] = (*newstate.T)[i];
      }
      nevLocal_ = blockSize_;
    }
    else {
      // get ritz vecs/vals
      Teuchos::SerialDenseMatrix<int,ScalarType> KK(blockSize_,blockSize_),
                                                 MM(blockSize_,blockSize_),
                                                  S(blockSize_,blockSize_);
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
#endif
        // project K
        MVT::MvTransMv(ONE,*X_,*KX_,KK);
        // project M
        MVT::MvTransMv(ONE,*X_,*MX_,MM);
        nevLocal_ = blockSize_;
      }

      // solve the projected problem
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerDS_ );
#endif
        Utils::directSolver(blockSize_, KK, Teuchos::rcpFromRef(MM), S, theta_, nevLocal_, 1);
        TEUCHOS_TEST_FOR_EXCEPTION(nevLocal_ != blockSize_,LOBPCGInitFailure,
                           "Anasazi::LOBPCG::initialize(): Initial Ritz analysis did not produce enough Ritz pairs to initialize algorithm.");
      }

      // We only have blockSize_ ritz pairs, ergo we do not need to select.
      // However, we still require them to be ordered correctly
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerSort_ );
#endif

        std::vector<int> order(blockSize_);
        // 
        // sort the first blockSize_ values in theta_
        sm_->sort(theta_, Teuchos::rcpFromRef(order), blockSize_);   // don't catch exception
        //
        // apply the same ordering to the primitive ritz vectors
        Utils::permuteVectors(order,S);
      }

      // update the solution, use R for storage
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
#endif
        // X <- X*S
        MVT::MvAddMv( ONE, *X_, ZERO, *X_, *R_ );        
        MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *X_ );
        // KX <- KX*S
        MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );        
        MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *KX_ );
        if (hasM_) {
          // MX <- MX*S
          MVT::MvAddMv( ONE, *MX_, ZERO, *MX_, *R_ );        
          MVT::MvTimesMatAddMv( ONE, *R_, S, ZERO, *MX_ );
        }
      }
    }

    //----------------------------------------
    // compute R
    //----------------------------------------
    if (newstate.R != Teuchos::null) {
      TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.R) != MVT::GetGlobalLength(*R_),
                          std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.R not correct." );
      TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) < blockSize_,
                          std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.R must have blockSize number of vectors." );
      MVT::SetBlock(*newstate.R,bsind,*R_);
    }
    else {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
#endif
      // form R <- KX - MX*T
      MVT::MvAddMv(ZERO,*KX_,ONE,*KX_,*R_);
      Teuchos::SerialDenseMatrix<int,ScalarType> T(blockSize_,blockSize_);
      for (int i=0; i<blockSize_; i++) T(i,i) = theta_[i];
      MVT::MvTimesMatAddMv(-ONE,*MX_,T,ONE,*R_);
    }

    // R has been updated; mark the norms as out-of-date
    Rnorms_current_ = false;
    R2norms_current_ = false;
  
    // put data in P,KP,MP: P is not used to set theta
    if (newstate.P != Teuchos::null) {
      TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.P) < blockSize_ ,
                          std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.P must have blockSize number of vectors." );
      TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.P) != MVT::GetGlobalLength(*P_),
                          std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.P not correct." );
      hasP_ = true;

      // set P_
      MVT::SetBlock(*newstate.P,bsind,*P_);

      // set/compute KP_
      if (newstate.KP != Teuchos::null) {
        TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.KP) < blockSize_,
                            std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.KP must have blockSize number of vectors." );
        TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.KP) != MVT::GetGlobalLength(*KP_),
                            std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.KP not correct." );
        MVT::SetBlock(*newstate.KP,bsind,*KP_);
      }
      else {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerOp_ );
#endif
        OPT::Apply(*Op_,*P_,*KP_);
        count_ApplyOp_ += blockSize_;
      }

      // set/compute MP_
      if (hasM_) {
        if (newstate.MP != Teuchos::null) {
          TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.MP) < blockSize_,
                              std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): newstate.MP must have blockSize number of vectors." );
          TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.MP) != MVT::GetGlobalLength(*MP_),
                              std::invalid_argument, "Anasazi::LOBPCG::initialize(newstate): vector length of newstate.MP not correct." );
          MVT::SetBlock(*newstate.MP,bsind,*MP_);
        }
        else {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
          Teuchos::TimeMonitor lcltimer( *timerMOp_ );
#endif
          OPT::Apply(*MOp_,*P_,*MP_);
          count_ApplyM_ += blockSize_;
        }
      }
    }
    else {
      hasP_ = false;
    }

    // finally, we are initialized
    initialized_ = true;

    if (om_->isVerbosity( Debug ) ) {
      // Check almost everything here
      CheckList chk;
      chk.checkX = true;
      chk.checkKX = true;
      chk.checkMX = true;
      chk.checkP = true;
      chk.checkKP = true;
      chk.checkMP = true;
      chk.checkR = true;
      chk.checkQ = true;
      om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
    }

  }

  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::initialize()
  {
    LOBPCGState<ScalarType,MV> empty;
    initialize(empty);
  }


  // Instruct the solver to use full orthogonalization
  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::setFullOrtho (bool fullOrtho)
  {
    if ( fullOrtho_ == true || initialized_ == false || fullOrtho == fullOrtho_ ) {
      // state is already orthogonalized or solver is not initialized
      fullOrtho_ = fullOrtho;
    }
    else {
      // solver is initialized, state is not fully orthogonalized, and user has requested full orthogonalization
      // ergo, we must throw away data in P
      fullOrtho_ = true;
      hasP_ = false;
    }

    // the user has called setFullOrtho, so the class has been instantiated
    // ergo, the data has already been allocated, i.e., setBlockSize() has been called
    // if it is already allocated, it should be the proper size
    if (fullOrtho_ && tmpmvec_ == Teuchos::null) {
      // allocated the workspace
      tmpmvec_ = MVT::Clone(*X_,blockSize_);
    }
    else if (fullOrtho_==false) {
      // free the workspace
      tmpmvec_ = Teuchos::null;
    }
  }


  // Perform LOBPCG iterations until the StatusTest tells us to stop.
  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::iterate () 
  {
    //
    // Allocate/initialize data structures
    //
    if (initialized_ == false) {
      initialize();
    }

    //
    // Miscellaneous definitions
    const int oneBlock    =   blockSize_;
    const int twoBlocks   = 2*blockSize_;
    const int threeBlocks = 3*blockSize_;

    std::vector<int> indblock1(blockSize_), indblock2(blockSize_), indblock3(blockSize_);
    for (int i=0; i<blockSize_; i++) {
      indblock1[i] = i;
      indblock2[i] = i +  blockSize_;
      indblock3[i] = i + 2*blockSize_;
    }

    //
    // Define dense projected/local matrices
    //   KK = Local stiffness matrix               (size: 3*blockSize_ x 3*blockSize_)
    //   MM = Local mass matrix                    (size: 3*blockSize_ x 3*blockSize_)
    //    S = Local eigenvectors                   (size: 3*blockSize_ x 3*blockSize_)
    Teuchos::SerialDenseMatrix<int,ScalarType> KK( threeBlocks, threeBlocks ), 
                                               MM( threeBlocks, threeBlocks ),
                                                S( threeBlocks, threeBlocks );

    while (tester_->checkStatus(this) != Passed) {

      // Print information on current status
      if (om_->isVerbosity(Debug)) {
        currentStatus( om_->stream(Debug) );
      }
      else if (om_->isVerbosity(IterationDetails)) {
        currentStatus( om_->stream(IterationDetails) );
      }

      // increment iteration counter
      iter_++;

      // Apply the preconditioner on the residuals: H <- Prec*R
      if (Prec_ != Teuchos::null) {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerPrec_ );
#endif
        OPT::Apply( *Prec_, *R_, *H_ );   // don't catch the exception
        count_ApplyPrec_ += blockSize_;
      }
      else {
        MVT::MvAddMv(ONE,*R_,ZERO,*R_,*H_);
      }

      // Apply the mass matrix on H
      if (hasM_) {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerMOp_ );
#endif
        OPT::Apply( *MOp_, *H_, *MH_);    // don't catch the exception
        count_ApplyM_ += blockSize_;
      }

      // orthogonalize H against the auxiliary vectors
      // optionally: orthogonalize H against X and P ([X P] is already orthonormal)
      Teuchos::Array<Teuchos::RCP<const MV> > Q;
      Q = auxVecs_;
      if (fullOrtho_) {
        // X and P are not contiguous, so there is not much point in putting them under 
        // a single multivector view
        Q.push_back(X_);
        if (hasP_) {
          Q.push_back(P_);
        }
      }
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerOrtho_ );
#endif
        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC = 
          Teuchos::tuple<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > >(Teuchos::null);
        int rank = orthman_->projectAndNormalizeMat(*H_,Q,dummyC,Teuchos::null,MH_);
        // our views are currently in place; it is safe to throw an exception
        TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,LOBPCGOrthoFailure,
                           "Anasazi::LOBPCG::iterate(): unable to compute orthonormal basis for H.");
      }

      if (om_->isVerbosity( Debug ) ) {
        CheckList chk;
        chk.checkH = true;
        chk.checkMH = true;
        om_->print( Debug, accuracyCheck(chk, ": after ortho H") );
      }
      else if (om_->isVerbosity( OrthoDetails ) ) {
        CheckList chk;
        chk.checkH = true;
        chk.checkMH = true;
        om_->print( OrthoDetails, accuracyCheck(chk,": after ortho H") );
      }

      // Apply the stiffness matrix to H
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerOp_ );
#endif
        OPT::Apply( *Op_, *H_, *KH_);   // don't catch the exception
        count_ApplyOp_ += blockSize_;
      }

      if (hasP_) {
        nevLocal_ = threeBlocks;
      }
      else {
        nevLocal_ = twoBlocks;
      }

      //
      // we need bases: [X H P] and [H P] (only need the latter if fullOrtho == false)
      // we need to perform the following operations:
      //    X' [KX KH KP]
      //    X' [MX MH MP]
      //    H' [KH KP]
      //    H' [MH MP]
      //    P' [KP]
      //    P' [MP]
      //    [X H P] CX
      //    [X H P] CP    if  fullOrtho
      //      [H P] CP    if !fullOrtho
      //
      // since M[X H P] is potentially the same memory as [X H P], and 
      // because we are not allowed to have overlapping non-const views of 
      // a multivector, we will now abandon our non-const views in favor of 
      // const views
      //
      X_ = Teuchos::null;
      KX_ = Teuchos::null;
      MX_ = Teuchos::null;
      H_ = Teuchos::null;
      KH_ = Teuchos::null;
      MH_ = Teuchos::null;
      P_ = Teuchos::null;
      KP_ = Teuchos::null;
      MP_ = Teuchos::null;
      Teuchos::RCP<const MV> cX, cH, cXHP, cHP, cK_XHP, cK_HP, cM_XHP, cM_HP, cP, cK_P, cM_P;
      {
        cX = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock1);
        cH = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock2);

        std::vector<int> indXHP(nevLocal_);
        for (int i=0; i<nevLocal_; i++) {
          indXHP[i] = i;
        }
        cXHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indXHP);
        cK_XHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indXHP);
        if (hasM_) {
          cM_XHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indXHP);
        }
        else {
          cM_XHP = cXHP;
        }

        std::vector<int> indHP(nevLocal_-blockSize_);
        for (int i=blockSize_; i<nevLocal_; i++) {
          indHP[i-blockSize_] = i;
        }
        cHP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indHP);
        cK_HP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indHP);
        if (hasM_) {
          cM_HP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indHP);
        }
        else {
          cM_HP = cHP;
        }

        if (nevLocal_ == threeBlocks) {
          cP = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(V_),indblock3);
          cK_P = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(KV_),indblock3);
          if (hasM_) {
            cM_P = MVT::CloneView(*Teuchos::rcp_implicit_cast<const MV>(MV_),indblock3);
          }
          else {
            cM_P = cP;
          }
        }
      }

      //
      //----------------------------------------
      // Form "local" mass and stiffness matrices
      //----------------------------------------
      {
        // We will form only the block upper triangular part of 
        // [X H P]' K [X H P]  and  [X H P]' M [X H P]
        // Get the necessary views into KK and MM:
        //      [--K1--]        [--M1--]
        // KK = [  -K2-]   MM = [  -M2-]
        //      [    K3]        [    M3]
        // 
        // It is okay to declare a zero-area view of a Teuchos::SerialDenseMatrix
        //
        Teuchos::SerialDenseMatrix<int,ScalarType> 
          K1(Teuchos::View,KK,blockSize_,nevLocal_             ,0*blockSize_,0*blockSize_),
          K2(Teuchos::View,KK,blockSize_,nevLocal_-1*blockSize_,1*blockSize_,1*blockSize_),
          K3(Teuchos::View,KK,blockSize_,nevLocal_-2*blockSize_,2*blockSize_,2*blockSize_),
          M1(Teuchos::View,MM,blockSize_,nevLocal_             ,0*blockSize_,0*blockSize_),
          M2(Teuchos::View,MM,blockSize_,nevLocal_-1*blockSize_,1*blockSize_,1*blockSize_),
          M3(Teuchos::View,MM,blockSize_,nevLocal_-2*blockSize_,2*blockSize_,2*blockSize_);
        {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
          Teuchos::TimeMonitor lcltimer( *timerLocalProj_ );
#endif
          MVT::MvTransMv( ONE, *cX, *cK_XHP, K1 );
          MVT::MvTransMv( ONE, *cX, *cM_XHP, M1 );
          MVT::MvTransMv( ONE, *cH, *cK_HP , K2 );
          MVT::MvTransMv( ONE, *cH, *cM_HP , M2 );
          if (nevLocal_ == threeBlocks) {
            MVT::MvTransMv( ONE, *cP, *cK_P, K3 );
            MVT::MvTransMv( ONE, *cP, *cM_P, M3 );
          }
        }
      }
      // below, we only need bases [X H P] and [H P] and friends
      // furthermore, we only need [H P] and friends if fullOrtho == false
      // clear the others now
      cX   = Teuchos::null;
      cH   = Teuchos::null;
      cP   = Teuchos::null;
      cK_P = Teuchos::null;
      cM_P = Teuchos::null;
      if (fullOrtho_ == true) {
        cHP = Teuchos::null;
        cK_HP = Teuchos::null;
        cM_HP = Teuchos::null;
      }

      //
      //---------------------------------------------------
      // Perform a spectral decomposition of (KK,MM)
      //---------------------------------------------------
      //
      // Get pointers to relevant part of KK, MM and S for Rayleigh-Ritz analysis
      Teuchos::SerialDenseMatrix<int,ScalarType> lclKK(Teuchos::View,KK,nevLocal_,nevLocal_), 
                                                 lclMM(Teuchos::View,MM,nevLocal_,nevLocal_),
                                                  lclS(Teuchos::View, S,nevLocal_,nevLocal_);
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerDS_ );
#endif
        int localSize = nevLocal_;
        Utils::directSolver(localSize, lclKK, Teuchos::rcpFromRef(lclMM), lclS, theta_, nevLocal_, 0);
        // localSize tells directSolver() how big KK,MM are
        // however, directSolver() may choose to use only the principle submatrices of KK,MM 
        // because of loss of MM-orthogonality in the projected eigenvectors
        // nevLocal_ tells us how much it used, telling us the effective localSize 
        // (i.e., how much of KK,MM used by directSolver)
        // we will not tolerate any indefiniteness, and will throw an exception if it was 
        // detected by directSolver
        //
        if (nevLocal_ != localSize) {
          // before throwing the exception, and thereby leaving iterate(), setup the views again
          // first, clear the const views
          cXHP   = Teuchos::null;
          cK_XHP = Teuchos::null;
          cM_XHP = Teuchos::null;
          cHP    = Teuchos::null;
          cK_HP  = Teuchos::null;
          cM_HP  = Teuchos::null;
          setupViews();
        }
        TEUCHOS_TEST_FOR_EXCEPTION(nevLocal_ != localSize, LOBPCGRitzFailure, 
            "Anasazi::LOBPCG::iterate(): indefiniteness detected in projected mass matrix." );
      }

      //
      //---------------------------------------------------
      // Sort the ritz values using the sort manager
      //---------------------------------------------------
      Teuchos::LAPACK<int,ScalarType> lapack;
      Teuchos::BLAS<int,ScalarType> blas;
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerSort_ );
#endif

        std::vector<int> order(nevLocal_);
        // 
        // Sort the first nevLocal_ values in theta_
        sm_->sort(theta_, Teuchos::rcpFromRef(order), nevLocal_);   // don't catch exception
        //
        // Sort the primitive ritz vectors
        Utils::permuteVectors(order,lclS);
      }

      //
      //----------------------------------------
      // Compute coefficients for X and P under [X H P]
      //----------------------------------------
      // Before computing X,P, optionally perform orthogonalization per Hetmaniuk,Lehoucq paper
      // CX will be the coefficients of [X,H,P] for new X, CP for new P
      // The paper suggests orthogonalizing CP against CX and orthonormalizing CP, w.r.t. MM
      // Here, we will also orthonormalize CX.
      // This is accomplished using the Cholesky factorization of [CX CP]^H lclMM [CX CP]
      Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > CX, CP;
      if (fullOrtho_) {
        // build orthonormal basis for (  0  ) that is MM orthogonal to ( S11 )
        //                             ( S21 )                          ( S21 )
        //                             ( S31 )                          ( S31 )
        // Do this using Cholesky factorization of ( S11  0  )
        //                                         ( S21 S21 )
        //                                         ( S31 S31 )
        //           ( S11  0  )
        // Build C = ( S21 S21 )
        //           ( S31 S31 )
        Teuchos::SerialDenseMatrix<int,ScalarType> C(nevLocal_,twoBlocks),
                                                MMC(nevLocal_,twoBlocks),
                                                CMMC(twoBlocks  ,twoBlocks);

        // first block of rows: ( S11 0 )
        for (int j=0; j<oneBlock; j++) {
          for (int i=0; i<oneBlock; i++) {
            // CX
            C(i,j) = lclS(i,j);
            // CP
            C(i,j+oneBlock) = ZERO;
          }
        }
        // second block of rows: (S21 S21)
        for (int j=0; j<oneBlock; j++) {
          for (int i=oneBlock; i<twoBlocks; i++) {
            // CX
            C(i,j)          = lclS(i,j);
            // CP
            C(i,j+oneBlock) = lclS(i,j);
          }
        }
        // third block of rows
        if (nevLocal_ == threeBlocks) {
          for (int j=0; j<oneBlock; j++) {
            for (int i=twoBlocks; i<threeBlocks; i++) {
              // CX
              C(i,j)          = lclS(i,j);
              // CP
              C(i,j+oneBlock) = lclS(i,j);
            }
          }
        }

        // compute MMC = lclMM*C
        {
          int teuchosret;
          teuchosret = MMC.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,C,ZERO);
          TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
              "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
          // compute CMMC = C^H*MMC == C^H*lclMM*C
          teuchosret = CMMC.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,C,MMC,ZERO);
          TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
              "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
        }

        // compute R (cholesky) of CMMC
        {
          int info;
          lapack.POTRF('U',twoBlocks,CMMC.values(),CMMC.stride(),&info);
          // our views ARE NOT currently in place; we must reestablish them before throwing an exception
          if (info != 0) {
            // Check symmetry of H'*K*H, this might be the first place a non-Hermitian operator will cause a problem.
            Teuchos::SerialDenseMatrix<int,ScalarType> K22(Teuchos::View,KK,blockSize_,blockSize_,1*blockSize_,1*blockSize_);
            Teuchos::SerialDenseMatrix<int,ScalarType> K22_trans( K22, Teuchos::CONJ_TRANS );
            K22_trans -= K22;
            MagnitudeType symNorm = K22_trans.normOne();
            MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));

            cXHP = Teuchos::null;
            cHP = Teuchos::null;
            cK_XHP = Teuchos::null;
            cK_HP = Teuchos::null;
            cM_XHP = Teuchos::null;
            cM_HP = Teuchos::null;
            setupViews();
           
            std::string errMsg = "Anasazi::LOBPCG::iterate(): Cholesky factorization failed during full orthogonalization.";
            if ( symNorm < tol )
            {
              TEUCHOS_TEST_FOR_EXCEPTION(info != 0, LOBPCGOrthoFailure, errMsg );
            }
            else
            {
              errMsg += " Check the operator A (or K), it may not be Hermitian!";
              TEUCHOS_TEST_FOR_EXCEPTION(info != 0, LOBPCGOrthoFailure, errMsg );
            }
          }
        }
        // compute C = C inv(R)
        blas.TRSM(Teuchos::RIGHT_SIDE,Teuchos::UPPER_TRI,Teuchos::NO_TRANS,Teuchos::NON_UNIT_DIAG,
                  nevLocal_,twoBlocks,ONE,CMMC.values(),CMMC.stride(),C.values(),C.stride());
        // put C(:,0:oneBlock-1) into CX
        CX = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,C,nevLocal_,oneBlock,0,0) );
        // put C(:,oneBlock:twoBlocks-1) into CP
        CP = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,C,nevLocal_,oneBlock,0,oneBlock) );

        // check the results
        if (om_->isVerbosity( Debug ) ) {
          Teuchos::SerialDenseMatrix<int,ScalarType> tmp1(nevLocal_,oneBlock),
                                                     tmp2(oneBlock,oneBlock);
          MagnitudeType tmp;
          int teuchosret;
          std::stringstream os;
          os.precision(2);
          os.setf(std::ios::scientific, std::ios::floatfield);

          os << " Checking Full Ortho: iteration " << iter_ << std::endl;

          // check CX^T MM CX == I
          // compute tmp1 = MM*CX
          teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CX,ZERO);
          TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
              "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
          // compute tmp2 = CX^H*tmp1 == CX^H*MM*CX
          teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CX,tmp1,ZERO);
          TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
              "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
          // subtrace tmp2 - I == CX^H * MM * CX - I
          for (int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
          tmp = tmp2.normFrobenius();          
          os << " >> Error in CX^H MM CX == I : " << tmp << std::endl;

          // check CP^T MM CP == I
          // compute tmp1 = MM*CP
          teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CP,ZERO);
          TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
              "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
          // compute tmp2 = CP^H*tmp1 == CP^H*MM*CP
          teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CP,tmp1,ZERO);
          TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,
              "Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
          // subtrace tmp2 - I == CP^H * MM * CP - I
          for (int i=0; i<oneBlock; i++) tmp2(i,i) -= ONE;
          tmp = tmp2.normFrobenius();          
          os << " >> Error in CP^H MM CP == I : " << tmp << std::endl;

          // check CX^T MM CP == 0
          // compute tmp1 = MM*CP
          teuchosret = tmp1.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,lclMM,*CP,ZERO);
          TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
          // compute tmp2 = CX^H*tmp1 == CX^H*MM*CP
          teuchosret = tmp2.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,*CX,tmp1,ZERO);
          TEUCHOS_TEST_FOR_EXCEPTION(teuchosret != 0,std::logic_error,"Anasazi::LOBPCG::iterate(): Logic error calling SerialDenseMatrix::multiply");
          // subtrace tmp2 == CX^H * MM * CP
          tmp = tmp2.normFrobenius();          
          os << " >> Error in CX^H MM CP == 0 : " << tmp << std::endl;

          os << std::endl;
          om_->print(Debug,os.str());
        }
      }
      else {
        //     [S11 ... ...]
        // S = [S21 ... ...]
        //     [S31 ... ...]
        //
        // CX = [S11]
        //      [S21]
        //      [S31]   ->  X = [X H P] CX
        //      
        // CP = [S21]   ->  P =   [H P] CP
        //      [S31]
        //
        CX = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,lclS,nevLocal_         ,oneBlock,0       ,0) );
        CP = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy,lclS,nevLocal_-oneBlock,oneBlock,oneBlock,0) );
      }

      //
      //----------------------------------------
      // Compute new X and new P
      //----------------------------------------
      // Note: Use R as a temporary work space and (if full ortho) tmpMV as well
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerLocalUpdate_ );
#endif

        // if full ortho, then CX and CP are dense
        // we multiply [X H P]*CX into tmpMV
        //             [X H P]*CP into R
        // then put V(:,firstblock) <- tmpMV
        //          V(:,thirdblock) <- R
        //
        // if no full ortho, then [H P]*CP doesn't reference first block (X) 
        // of V, so that we can modify it before computing P
        // so we multiply [X H P]*CX into R
        //                V(:,firstblock) <- R
        //       multiply [H P]*CP into R
        //                V(:,thirdblock) <- R
        //
        // mutatis mutandis for K[XP] and M[XP]
        //
        // use SetBlock to do the assignments into V_
        //
        // in either case, views are only allowed to be overlapping
        // if they are const, and it should be assume that SetBlock
        // creates a view of the associated part
        //
        // we have from above const-pointers to [KM]XHP, [KM]HP and (if hasP) [KM]P
        //
        if (fullOrtho_) {
          // X,P
          MVT::MvTimesMatAddMv(ONE,*cXHP,*CX,ZERO,*tmpmvec_);
          MVT::MvTimesMatAddMv(ONE,*cXHP,*CP,ZERO,*R_);
          cXHP = Teuchos::null;
          MVT::SetBlock(*tmpmvec_,indblock1,*V_);
          MVT::SetBlock(*R_      ,indblock3,*V_);
          // KX,KP
          MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CX,ZERO,*tmpmvec_);
          MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CP,ZERO,*R_);
          cK_XHP = Teuchos::null;
          MVT::SetBlock(*tmpmvec_,indblock1,*KV_);
          MVT::SetBlock(*R_      ,indblock3,*KV_);
          // MX,MP
          if (hasM_) {
            MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CX,ZERO,*tmpmvec_);
            MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CP,ZERO,*R_);
            cM_XHP = Teuchos::null;
            MVT::SetBlock(*tmpmvec_,indblock1,*MV_);
            MVT::SetBlock(*R_      ,indblock3,*MV_);
          }
          else {
            cM_XHP = Teuchos::null;
          }
        }
        else {
          // X,P
          MVT::MvTimesMatAddMv(ONE,*cXHP,*CX,ZERO,*R_);
          cXHP = Teuchos::null;
          MVT::SetBlock(*R_,indblock1,*V_);
          MVT::MvTimesMatAddMv(ONE,*cHP,*CP,ZERO,*R_);
          cHP = Teuchos::null;
          MVT::SetBlock(*R_,indblock3,*V_);
          // KX,KP
          MVT::MvTimesMatAddMv(ONE,*cK_XHP,*CX,ZERO,*R_);
          cK_XHP = Teuchos::null;
          MVT::SetBlock(*R_,indblock1,*KV_);
          MVT::MvTimesMatAddMv(ONE,*cK_HP,*CP,ZERO,*R_);
          cK_HP = Teuchos::null;
          MVT::SetBlock(*R_,indblock3,*KV_);
          // MX,MP
          if (hasM_) {
            MVT::MvTimesMatAddMv(ONE,*cM_XHP,*CX,ZERO,*R_);
            cM_XHP = Teuchos::null;
            MVT::SetBlock(*R_,indblock1,*MV_);
            MVT::MvTimesMatAddMv(ONE,*cM_HP,*CP,ZERO,*R_);
            cM_HP = Teuchos::null;
            MVT::SetBlock(*R_,indblock3,*MV_);
          }
          else {
            cM_XHP = Teuchos::null;
            cM_HP = Teuchos::null;
          }
        }
      } // end timing block
      // done with coefficient matrices
      CX = Teuchos::null;
      CP = Teuchos::null;

      //
      // we now have a P direction
      hasP_ = true;

      // debugging check: all of our const views should have been cleared by now
      // if not, we have a logic error above
      TEUCHOS_TEST_FOR_EXCEPTION(   cXHP != Teuchos::null || cK_XHP != Teuchos::null || cM_XHP != Teuchos::null
                          || cHP != Teuchos::null ||  cK_HP != Teuchos::null || cM_HP  != Teuchos::null
                          ||  cP != Teuchos::null ||   cK_P != Teuchos::null || cM_P   != Teuchos::null,
                          std::logic_error,
                          "Anasazi::BlockKrylovSchur::iterate(): const views were not all cleared! Something went wrong!" );

      //
      // recreate our const MV views of X,H,P and friends
      setupViews();

      //
      // Compute the new residuals, explicitly
      {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
        Teuchos::TimeMonitor lcltimer( *timerCompRes_ );
#endif
        MVT::MvAddMv( ONE, *KX_, ZERO, *KX_, *R_ );
        Teuchos::SerialDenseMatrix<int,ScalarType> T( blockSize_, blockSize_ );
        for (int i = 0; i < blockSize_; i++) {
          T(i,i) = theta_[i];
        }
        MVT::MvTimesMatAddMv( -ONE, *MX_, T, ONE, *R_ );
      }

      // R has been updated; mark the norms as out-of-date
      Rnorms_current_ = false;
      R2norms_current_ = false;

      // When required, monitor some orthogonalities
      if (om_->isVerbosity( Debug ) ) {
        // Check almost everything here
        CheckList chk;
        chk.checkX = true;
        chk.checkKX = true;
        chk.checkMX = true;
        chk.checkP = true;
        chk.checkKP = true;
        chk.checkMP = true;
        chk.checkR = true;
        om_->print( Debug, accuracyCheck(chk, ": after local update") );
      }
      else if (om_->isVerbosity( OrthoDetails )) {
        CheckList chk;
        chk.checkX = true;
        chk.checkP = true;
        chk.checkR = true;
        om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
      }
    } // end while (statusTest == false)
  }


  // compute/return residual M-norms
  template <class ScalarType, class MV, class OP>
  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
  LOBPCG<ScalarType,MV,OP>::getResNorms() {
    if (Rnorms_current_ == false) {
      // Update the residual norms
      orthman_->norm(*R_,Rnorms_);
      Rnorms_current_ = true;
    }
    return Rnorms_;
  }

  
  // compute/return residual 2-norms
  template <class ScalarType, class MV, class OP>
  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
  LOBPCG<ScalarType,MV,OP>::getRes2Norms() {
    if (R2norms_current_ == false) {
      // Update the residual 2-norms 
      MVT::MvNorm(*R_,R2norms_);
      R2norms_current_ = true;
    }
    return R2norms_;
  }




  // Check accuracy, orthogonality, and other debugging stuff
  // 
  // bools specify which tests we want to run (instead of running more than we actually care about)
  //
  // we don't bother checking the following because they are computed explicitly:
  //    H == Prec*R
  //   KH == K*H
  //
  // 
  // checkX : X orthonormal
  //          orthogonal to auxvecs
  // checkMX: check MX == M*X
  // checkKX: check KX == K*X
  // checkP : if fullortho P orthonormal and orthogonal to X
  //          orthogonal to auxvecs
  // checkMP: check MP == M*P
  // checkKP: check KP == K*P
  // checkH : if fullortho H orthonormal and orthogonal to X and P
  //          orthogonal to auxvecs
  // checkMH: check MH == M*H
  // checkR : check R orthogonal to X
  // checkQ : check that auxiliary vectors are actually orthonormal
  //
  // TODO: 
  //  add checkTheta 
  //
  template <class ScalarType, class MV, class OP>
  std::string LOBPCG<ScalarType,MV,OP>::accuracyCheck( const CheckList &chk, const std::string &where ) const 
  {
    using std::endl;

    std::stringstream os;
    os.precision(2);
    os.setf(std::ios::scientific, std::ios::floatfield);
    MagnitudeType tmp;

    os << " Debugging checks: iteration " << iter_ << where << endl;

    // X and friends
    if (chk.checkX && initialized_) {
      tmp = orthman_->orthonormError(*X_);
      os << " >> Error in X^H M X == I : " << tmp << endl;
      for (Array_size_type i=0; i<auxVecs_.size(); i++) {
        tmp = orthman_->orthogError(*X_,*auxVecs_[i]);
        os << " >> Error in X^H M Q[" << i << "] == 0 : " << tmp << endl;
      }
    }
    if (chk.checkMX && hasM_ && initialized_) {
      tmp = Utils::errorEquality(*X_, *MX_, MOp_);
      os << " >> Error in MX == M*X    : " << tmp << endl;
    }
    if (chk.checkKX && initialized_) {
      tmp = Utils::errorEquality(*X_, *KX_, Op_);
      os << " >> Error in KX == K*X    : " << tmp << endl;
    }

    // P and friends
    if (chk.checkP && hasP_ && initialized_) {
      if (fullOrtho_) {
        tmp = orthman_->orthonormError(*P_);
        os << " >> Error in P^H M P == I : " << tmp << endl;
        tmp = orthman_->orthogError(*P_,*X_);
        os << " >> Error in P^H M X == 0 : " << tmp << endl;
      }
      for (Array_size_type i=0; i<auxVecs_.size(); i++) {
        tmp = orthman_->orthogError(*P_,*auxVecs_[i]);
        os << " >> Error in P^H M Q[" << i << "] == 0 : " << tmp << endl;
      }
    }
    if (chk.checkMP && hasM_ && hasP_ && initialized_) {
      tmp = Utils::errorEquality(*P_, *MP_, MOp_);
      os << " >> Error in MP == M*P    : " << tmp << endl;
    }
    if (chk.checkKP && hasP_ && initialized_) {
      tmp = Utils::errorEquality(*P_, *KP_, Op_);
      os << " >> Error in KP == K*P    : " << tmp << endl;
    }

    // H and friends
    if (chk.checkH && initialized_) {
      if (fullOrtho_) {
        tmp = orthman_->orthonormError(*H_);
        os << " >> Error in H^H M H == I : " << tmp << endl;
        tmp = orthman_->orthogError(*H_,*X_);
        os << " >> Error in H^H M X == 0 : " << tmp << endl;
        if (hasP_) {
          tmp = orthman_->orthogError(*H_,*P_);
          os << " >> Error in H^H M P == 0 : " << tmp << endl;
        }
      }
      for (Array_size_type i=0; i<auxVecs_.size(); i++) {
        tmp = orthman_->orthogError(*H_,*auxVecs_[i]);
        os << " >> Error in H^H M Q[" << i << "] == 0 : " << tmp << endl;
      }
    }
    if (chk.checkMH && hasM_ && initialized_) {
      tmp = Utils::errorEquality(*H_, *MH_, MOp_);
      os << " >> Error in MH == M*H    : " << tmp << endl;
    }

    // R: this is not M-orthogonality, but standard euclidean orthogonality
    if (chk.checkR && initialized_) {
      Teuchos::SerialDenseMatrix<int,ScalarType> xTx(blockSize_,blockSize_);
      MVT::MvTransMv(ONE,*X_,*R_,xTx);
      tmp = xTx.normFrobenius();
      MVT::MvTransMv(ONE,*R_,*R_,xTx);
      double normR = xTx.normFrobenius();
      os << " >> RelError in X^H R == 0: " << tmp/normR << endl;
    }

    // Q
    if (chk.checkQ) {
      for (Array_size_type i=0; i<auxVecs_.size(); i++) {
        tmp = orthman_->orthonormError(*auxVecs_[i]);
        os << " >> Error in Q[" << i << "]^H M Q[" << i << "] == I : " << tmp << endl;
        for (Array_size_type j=i+1; j<auxVecs_.size(); j++) {
          tmp = orthman_->orthogError(*auxVecs_[i],*auxVecs_[j]);
          os << " >> Error in Q[" << i << "]^H M Q[" << j << "] == 0 : " << tmp << endl;
        }
      }
    }

    os << endl;

    return os.str();
  }


  // Print the current status of the solver
  template <class ScalarType, class MV, class OP>
  void 
  LOBPCG<ScalarType,MV,OP>::currentStatus(std::ostream &os) 
  {
    using std::endl;

    os.setf(std::ios::scientific, std::ios::floatfield);  
    os.precision(6);
    os <<endl;
    os <<"================================================================================" << endl;
    os << endl;
    os <<"                              LOBPCG Solver Status" << endl;
    os << endl;
    os <<"The solver is "<<(initialized_ ? "initialized." : "not initialized.") << endl;
    os <<"The number of iterations performed is " << iter_       << endl;
    os <<"The current block size is             " << blockSize_  << endl;
    os <<"The number of auxiliary vectors is    " << numAuxVecs_ << endl;
    os <<"The number of operations Op*x   is " << count_ApplyOp_   << endl;
    os <<"The number of operations M*x    is " << count_ApplyM_    << endl;
    os <<"The number of operations Prec*x is " << count_ApplyPrec_ << endl;

    os.setf(std::ios_base::right, std::ios_base::adjustfield);

    if (initialized_) {
      os << endl;
      os <<"CURRENT EIGENVALUE ESTIMATES             "<<endl;
      os << std::setw(20) << "Eigenvalue" 
         << std::setw(20) << "Residual(M)"
         << std::setw(20) << "Residual(2)"
         << endl;
      os <<"--------------------------------------------------------------------------------"<<endl;
      for (int i=0; i<blockSize_; i++) {
        os << std::setw(20) << theta_[i];
        if (Rnorms_current_) os << std::setw(20) << Rnorms_[i];
        else os << std::setw(20) << "not current";
        if (R2norms_current_) os << std::setw(20) << R2norms_[i];
        else os << std::setw(20) << "not current";
        os << endl;
      }
    }
    os <<"================================================================================" << endl;
    os << endl;
  }

  // are we initialized or not?
  template <class ScalarType, class MV, class OP>
  bool LOBPCG<ScalarType,MV,OP>::isInitialized() const { 
    return initialized_; 
  }


  // is P valid or not?
  template <class ScalarType, class MV, class OP>
  bool LOBPCG<ScalarType,MV,OP>::hasP() {
    return hasP_;
  }
  
  // is full orthogonalization enabled or not?
  template <class ScalarType, class MV, class OP>
  bool LOBPCG<ScalarType,MV,OP>::getFullOrtho() const { 
    return(fullOrtho_); 
  }

  
  // return the current auxilliary vectors
  template <class ScalarType, class MV, class OP>
  Teuchos::Array<Teuchos::RCP<const MV> > LOBPCG<ScalarType,MV,OP>::getAuxVecs() const {
    return auxVecs_;
  }

  // return the current block size
  template <class ScalarType, class MV, class OP>
  int LOBPCG<ScalarType,MV,OP>::getBlockSize() const {
    return(blockSize_); 
  }

  // return the current eigenproblem
  template <class ScalarType, class MV, class OP>
  const Eigenproblem<ScalarType,MV,OP>& LOBPCG<ScalarType,MV,OP>::getProblem() const { 
    return(*problem_); 
  }

  
  // return the max subspace dimension
  template <class ScalarType, class MV, class OP>
  int LOBPCG<ScalarType,MV,OP>::getMaxSubspaceDim() const {
    return 3*blockSize_;
  }

  // return the current subspace dimension
  template <class ScalarType, class MV, class OP>
  int LOBPCG<ScalarType,MV,OP>::getCurSubspaceDim() const {
    if (!initialized_) return 0;
    return nevLocal_;
  }

  
  // return the current ritz residual norms
  template <class ScalarType, class MV, class OP>
  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> 
  LOBPCG<ScalarType,MV,OP>::getRitzRes2Norms() 
  {
    return this->getRes2Norms();
  }

  
  // return the current compression indices
  template <class ScalarType, class MV, class OP>
  std::vector<int> LOBPCG<ScalarType,MV,OP>::getRitzIndex() {
    std::vector<int> ret(nevLocal_,0);
    return ret;
  }

  
  // return the current ritz values
  template <class ScalarType, class MV, class OP>
  std::vector<Value<ScalarType> > LOBPCG<ScalarType,MV,OP>::getRitzValues() { 
    std::vector<Value<ScalarType> > ret(nevLocal_);
    for (int i=0; i<nevLocal_; i++) {
      ret[i].realpart = theta_[i];
      ret[i].imagpart = ZERO;
    }
    return ret;
  }

  // Set a new StatusTest for the solver.
  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::setStatusTest(Teuchos::RCP<StatusTest<ScalarType,MV,OP> > test) {
    TEUCHOS_TEST_FOR_EXCEPTION(test == Teuchos::null,std::invalid_argument,
        "Anasazi::LOBPCG::setStatusTest() was passed a null StatusTest.");
    tester_ = test;
  }

  // Get the current StatusTest used by the solver.
  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > LOBPCG<ScalarType,MV,OP>::getStatusTest() const {
    return tester_;
  }
  
  // return the current ritz vectors
  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<const MV> LOBPCG<ScalarType,MV,OP>::getRitzVectors() {
    return X_;
  }

  
  // reset the iteration counter
  template <class ScalarType, class MV, class OP>
  void LOBPCG<ScalarType,MV,OP>::resetNumIters() {
    iter_=0; 
  }

  
  // return the number of iterations
  template <class ScalarType, class MV, class OP>
  int LOBPCG<ScalarType,MV,OP>::getNumIters() const { 
    return(iter_); 
  }

  
  // return the state
  template <class ScalarType, class MV, class OP>
  LOBPCGState<ScalarType,MV> LOBPCG<ScalarType,MV,OP>::getState() const {
    LOBPCGState<ScalarType,MV> state;
    state.V = V_;
    state.KV = KV_;
    state.X = X_;
    state.KX = KX_;
    state.P = P_;
    state.KP = KP_;
    state.H = H_;
    state.KH = KH_;
    state.R = R_;
    state.T = Teuchos::rcp(new std::vector<MagnitudeType>(theta_));
    if (hasM_) {
      state.MV = MV_;
      state.MX = MX_;
      state.MP = MP_;
      state.MH = MH_;
    }
    else {
      state.MX = Teuchos::null;
      state.MP = Teuchos::null;
      state.MH = Teuchos::null;
    }
    return state;
  }

} // end Anasazi namespace

#endif // ANASAZI_LOBPCG_HPP