AnasaziBasicEigenproblem.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                 Anasazi: Block Eigensolvers Package
//                 Copyright 2004 Sandia Corporation
//
// Under 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 ANASAZI_BASIC_EIGENPROBLEM_H
#define ANASAZI_BASIC_EIGENPROBLEM_H

#include "AnasaziEigenproblem.hpp"
#include "AnasaziMultiVecTraits.hpp"
#include "AnasaziOperatorTraits.hpp"

namespace Anasazi {
  
  template<class ScalarType, class MV, class OP>
  class BasicEigenproblem : public virtual Eigenproblem<ScalarType, MV, OP> {
    
  public:
    

    
    BasicEigenproblem();
    
    BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<MV>& InitVec );
    
    BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<const OP>& B, const Teuchos::RCP<MV>& InitVec );
    
    BasicEigenproblem( const BasicEigenproblem<ScalarType, MV, OP>& Problem );
    
    virtual ~BasicEigenproblem() {};
    

    
    void setOperator( const Teuchos::RCP<const OP>& Op ) { _Op = Op; _isSet=false; };
    
    void setA( const Teuchos::RCP<const OP>& A ) { _AOp = A; _isSet=false; };
    
    void setM( const Teuchos::RCP<const OP>& M ) { _MOp = M; _isSet=false; };
    
    void setPrec( const Teuchos::RCP<const OP>& Prec ) { _Prec = Prec; _isSet=false; };
    
    void setInitVec( const Teuchos::RCP<MV>& InitVec ) { _InitVec = InitVec; _isSet=false; };
    
    void setAuxVecs( const Teuchos::RCP<const MV>& AuxVecs ) { _AuxVecs = AuxVecs; _isSet=false; };

    void setNEV( int nev ){ _nev = nev; _isSet=false; };


    void setHermitian( bool isSym ){ _isSym = isSym; _isSet=false; };

    bool setProblem();

    void setSolution(const Eigensolution<ScalarType,MV> &sol) {_sol = sol;}

    

    
    Teuchos::RCP<const OP> getOperator() const { return( _Op ); };
    
    Teuchos::RCP<const OP> getA() const { return( _AOp ); };
    
    Teuchos::RCP<const OP> getM() const { return( _MOp ); };
    
    Teuchos::RCP<const OP> getPrec() const { return( _Prec ); };
    
    Teuchos::RCP<const MV> getInitVec() const { return( _InitVec ); };
    
    Teuchos::RCP<const MV> getAuxVecs() const { return( _AuxVecs ); };
    
    int getNEV() const { return( _nev ); }

    bool isHermitian() const { return( _isSym ); }
    
    bool isProblemSet() const { return( _isSet ); }

    const Eigensolution<ScalarType,MV> & getSolution() const { return(_sol); }

    


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


  protected:
    
    Teuchos::RCP<const OP> _AOp;

    Teuchos::RCP<const OP> _MOp; 

    Teuchos::RCP<const OP> _Op;

    Teuchos::RCP<const OP> _Prec;

    Teuchos::RCP<MV> _InitVec;

    Teuchos::RCP<const MV> _AuxVecs;

    int _nev;


    bool _isSym;

    bool _isSet;

    typedef MultiVecTraits<ScalarType,MV> MVT;
    typedef OperatorTraits<ScalarType,MV,OP> OPT;

    Eigensolution<ScalarType,MV> _sol;
  };


  //=============================================================================
  //     Implementations (Constructors / Destructors)
  //=============================================================================
  template <class ScalarType, class MV, class OP>
  BasicEigenproblem<ScalarType, MV, OP>::BasicEigenproblem() : 
    _nev(0), 
    _isSym(false),
    _isSet(false)
  {
  }


  //=============================================================================
  template <class ScalarType, class MV, class OP>
  BasicEigenproblem<ScalarType, MV, OP>::BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<MV>& InitVec ) :    
    _Op(Op), 
    _InitVec(InitVec), 
    _nev(0), 
    _isSym(false),
    _isSet(false)
  {
  }


  //=============================================================================
  template <class ScalarType, class MV, class OP>
  BasicEigenproblem<ScalarType, MV, OP>::BasicEigenproblem( const Teuchos::RCP<const OP>& Op, const Teuchos::RCP<const OP>& M,
                                                            const Teuchos::RCP<MV>& InitVec ) :
    _MOp(M), 
    _Op(Op), 
    _InitVec(InitVec), 
    _nev(0), 
    _isSym(false),
    _isSet(false)
  {
  }


  //=============================================================================
  template <class ScalarType, class MV, class OP>
  BasicEigenproblem<ScalarType, MV, OP>::BasicEigenproblem( const BasicEigenproblem<ScalarType,MV,OP>& Problem ) :
    _AOp(Problem._AOp), 
    _MOp(Problem._MOp), 
    _Op(Problem._Op), 
    _Prec(Problem._Prec), 
    _InitVec(Problem._InitVec), 
    _nev(Problem._nev), 
    _isSym(Problem._isSym),
    _isSet(Problem._isSet),
    _sol(Problem._sol)
  {
  }


  //=============================================================================
  //    SetProblem (sanity check method)
  //=============================================================================
  template <class ScalarType, class MV, class OP>
  bool BasicEigenproblem<ScalarType, MV, OP>::setProblem() 
  {
    //----------------------------------------------------------------
    // Sanity Checks
    //----------------------------------------------------------------
    // If there is no operator, then we can't proceed.
    if ( !_AOp.get() && !_Op.get() ) { return false; }
    
    // If there is no initial vector, then we don't have anything to clone workspace from.
    if ( !_InitVec.get() ) { return false; }
    
    // If we don't need any eigenvalues, we don't need to continue.
    if (_nev == 0) { return false; }
    
    // If there is an A, but no operator, we can set them equal.
    if (_AOp.get() && !_Op.get()) { _Op = _AOp; }

    // Clear the storage from any previous call to setSolution()
    Eigensolution<ScalarType,MV> emptysol;
    _sol = emptysol;
    
    // mark the problem as set and return no-error
    _isSet=true;
    return true;
  }

  //=============================================================================
  //    Computes the residual vector
  //=============================================================================
  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<const MV> BasicEigenproblem<ScalarType, MV, OP>::computeCurrResVec() const
  {
    using Teuchos::RCP;

    TEUCHOS_TEST_FOR_EXCEPTION(!isHermitian(), std::invalid_argument,
        "BasicEigenproblem::computeCurrResVec: This method is not currently "
        "implemented for non-Hermitian problems.  Sorry for any inconvenience.");

    const Eigensolution<ScalarType,MV> sol = getSolution();

    if(sol.numVecs <= 0)
      return Teuchos::null;

    // Copy the eigenvalues
    RCP<MV> X = sol.Evecs;
    std::vector<ScalarType> Lambda(sol.numVecs);
    for(int i = 0; i < sol.numVecs; i++)
      Lambda[i] = sol.Evals[i].realpart;

    // Compute AX
    RCP<MV> AX = MVT::Clone(*X,sol.numVecs);
    if(_AOp != Teuchos::null)
      OPT::Apply(*_AOp,*X,*AX);
    else
      OPT::Apply(*_Op,*X,*AX);

    // Compute MX if necessary
    RCP<MV> MX;
    if(_MOp != Teuchos::null)
    {
      MX = MVT::Clone(*X,sol.numVecs);
      OPT::Apply(*_MOp,*X,*MX);
    }
    else
    {
      MX = MVT::CloneCopy(*X);
    }

    // Compute R = AX - MX \Lambda
    RCP<MV> R = MVT::Clone(*X,sol.numVecs);
    MVT::MvScale(*MX,Lambda);
    MVT::MvAddMv(1.0,*AX,-1.0,*MX,*R);

    return R;
  }
  
} // end Anasazi namespace
#endif

// end AnasaziBasicEigenproblem.hpp