AnasaziTraceMinDavidsonSolMgr.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                 Anasazi: Block Eigensolvers Package
//                 Copyright (2004) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER

#ifndef ANASAZI_TRACEMIN_DAVIDSON_SOLMGR_HPP
#define ANASAZI_TRACEMIN_DAVIDSON_SOLMGR_HPP

#include "AnasaziConfigDefs.hpp"
#include "AnasaziTraceMinBaseSolMgr.hpp"
#include "AnasaziTraceMinDavidson.hpp"

namespace Anasazi {
namespace Experimental {

template<class ScalarType, class MV, class OP>
class TraceMinDavidsonSolMgr : public TraceMinBaseSolMgr<ScalarType,MV,OP> {

  private:
    typedef MultiVecTraits<ScalarType,MV> MVT;
    typedef MultiVecTraitsExt<ScalarType,MV> MVText;
    typedef OperatorTraits<ScalarType,MV,OP> OPT;
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
    typedef Teuchos::ScalarTraits<MagnitudeType> MT;
    
  public:



  TraceMinDavidsonSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
                          Teuchos::ParameterList &pl );

  private:
    int maxRestarts_;

    // Returns true if the subspace is full
    bool needToRestart(const Teuchos::RCP< TraceMinBase<ScalarType,MV,OP> > solver)
    { 
      return (solver->getCurSubspaceDim() == solver->getMaxSubspaceDim()); 
    };

    // Performs a restart by reinitializing TraceMinDavidson with the most significant part of the basis
    bool performRestart(int &numRestarts, Teuchos::RCP< TraceMinBase<ScalarType,MV,OP> > solver);

    // Returns a TraceMinDavidson solver
    Teuchos::RCP< TraceMinBase<ScalarType,MV,OP> > createSolver( 
            const Teuchos::RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
            const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >      &outputtest,
            const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
            Teuchos::ParameterList &plist)
    {
      return Teuchos::rcp( new TraceMinDavidson<ScalarType,MV,OP>(this->problem_,sorter,this->printer_,outputtest,ortho,plist) );
    };
};


// Basic constructor for TraceMinDavidsonSolMgr
template<class ScalarType, class MV, class OP>
TraceMinDavidsonSolMgr<ScalarType,MV,OP>::TraceMinDavidsonSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem, Teuchos::ParameterList &pl ) :
  TraceMinBaseSolMgr<ScalarType,MV,OP>(problem,pl)
{
  // TODO: Come back tot these exceptions and make the descriptions better.
  maxRestarts_ = pl.get("Maximum Restarts", 50);
  TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts_ <= 0, std::invalid_argument,
         "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Maximum Restarts\" must be strictly positive.");

  // block size: default is 1
  this->blockSize_ = pl.get("Block Size", 1);
  TEUCHOS_TEST_FOR_EXCEPTION(this->blockSize_ <= 0, std::invalid_argument,
         "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Block Size\" must be strictly positive.");

  this->numRestartBlocks_ = (int)std::ceil(this->problem_->getNEV() / this->blockSize_);
  this->numRestartBlocks_ = pl.get("Num Restart Blocks", this->numRestartBlocks_);
  TEUCHOS_TEST_FOR_EXCEPTION(this->numRestartBlocks_ <= 0, std::invalid_argument,
         "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Num Restart Blocks\" must be strictly positive.");

  this->numBlocks_ = pl.get("Num Blocks", 3*this->numRestartBlocks_);
  TEUCHOS_TEST_FOR_EXCEPTION(this->numBlocks_ <= 1, std::invalid_argument,
         "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Num Blocks\" must be greater than 1.  If you only wish to use one block, please use TraceMinSolMgr instead of TraceMinDavidsonSolMgr.");

  TEUCHOS_TEST_FOR_EXCEPTION(this->numRestartBlocks_ >= this->numBlocks_, std::invalid_argument,
         "Anasazi::TraceMinDavidsonSolMgr::constructor(): \"Num Blocks\" must be strictly greater than \"Num Restart Blocks\".");

  std::stringstream ss;
  ss << "Anasazi::TraceMinDavidsonSolMgr::constructor(): Potentially impossible orthogonality requests. Reduce basis size (" << static_cast<ptrdiff_t>(this->numBlocks_)*this->blockSize_ << ") or locking size (" << this->maxLocked_ << ") because " << static_cast<ptrdiff_t>(this->numBlocks_) << "*" << this->blockSize_ << " + " << this->maxLocked_ << " > " << MVText::GetGlobalLength(*this->problem_->getInitVec()) << ".";
  TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(this->numBlocks_)*this->blockSize_ + this->maxLocked_ > MVText::GetGlobalLength(*this->problem_->getInitVec()),
         std::invalid_argument, ss.str());

  TEUCHOS_TEST_FOR_EXCEPTION(this->maxLocked_ + this->blockSize_ < this->problem_->getNEV(), std::invalid_argument,
         "Anasazi::TraceMinDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
}


// Performs a restart by reinitializing TraceMinDavidson with the most significant part of the basis
template <class ScalarType, class MV, class OP>
bool TraceMinDavidsonSolMgr<ScalarType,MV,OP>::performRestart(int &numRestarts, Teuchos::RCP< TraceMinBase<ScalarType,MV,OP> > solver)
{
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
  Teuchos::TimeMonitor restimer(*this->_timerRestarting);
#endif

  if ( numRestarts >= maxRestarts_ ) {
    return false; // break from while(1){tm_solver->iterate()}
  }
  numRestarts++;

  this->printer_->stream(IterationDetails) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;

  TraceMinBaseState<ScalarType,MV> oldstate = solver->getState();
  TraceMinBaseState<ScalarType,MV> newstate;
  int newdim = this->numRestartBlocks_*this->blockSize_;
  std::vector<int> indToCopy(newdim);
  for(int i=0; i<newdim; i++) indToCopy[i] = i;

  // Copy the relevant parts of the old state to the new one
  // This may involve computing parts of X
  this->copyPartOfState (oldstate, newstate, indToCopy);

  // send the new state to the solver
  newstate.NEV = oldstate.NEV;
  solver->initialize(newstate); 

  return true;
}


}} // end Anasazi namespace

#endif /* ANASAZI_TraceMinDavidson_SOLMGR_HPP */