AnasaziTraceMin.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
//
//                 Anasazi: Block Eigensolvers Package
//

#ifndef ANASAZI_TRACEMIN_HPP
#define ANASAZI_TRACEMIN_HPP

#include "AnasaziTypes.hpp"
#include "AnasaziBasicSort.hpp"
#include "AnasaziTraceMinBase.hpp"

#ifdef HAVE_ANASAZI_EPETRA
  #include "Epetra_Operator.h"
#endif

#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_SerialDenseSolver.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_TimeMonitor.hpp"

// TODO: TraceMin performs some unnecessary computations upon restarting.  Fix it!

namespace Anasazi {
namespace Experimental {
  template <class ScalarType, class MV, class OP>
  class TraceMin : public TraceMinBase<ScalarType,MV,OP> {
  public:


    TraceMin( 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
                 );

  private:
    //
    // 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;

    // TraceMin specific methods
    void addToBasis(const Teuchos::RCP<const MV> Delta);

    void harmonicAddToBasis(const Teuchos::RCP<const MV> Delta);
  };

  //
  // Implementations
  //


  // Constructor
  template <class ScalarType, class MV, class OP>
  TraceMin<ScalarType,MV,OP>::TraceMin(
        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
        ) :
    TraceMinBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params),
    ONE(Teuchos::ScalarTraits<MagnitudeType>::one()),
    ZERO(Teuchos::ScalarTraits<MagnitudeType>::zero()),
    NANVAL(Teuchos::ScalarTraits<MagnitudeType>::nan())
  {
  }


  template <class ScalarType, class MV, class OP>
  void TraceMin<ScalarType,MV,OP>::addToBasis(const Teuchos::RCP<const MV> Delta)
  {
    MVT::MvAddMv(ONE,*this->X_,-ONE,*Delta,*this->V_);

    if(this->hasM_)
    {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
#endif
      this->count_ApplyM_+= this->blockSize_;

      OPT::Apply(*this->MOp_,*this->V_,*this->MV_);
    }

    int rank;
    {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
#endif

      if(this->numAuxVecs_ > 0)
      {
        rank = this->orthman_->projectAndNormalizeMat(*this->V_,this->auxVecs_,
               Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)),
               Teuchos::null,this->MV_,this->MauxVecs_);
      }
      else
      {
        rank = this->orthman_->normalizeMat(*this->V_,Teuchos::null,this->MV_);
      }
    }

    // FIXME (mfh 07 Oct 2014) This variable is currently unused, but
    // it would make sense to use it to check whether the block is
    // rank deficient.
    (void) rank;

    if(this->Op_ != Teuchos::null)
    {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
#endif
      this->count_ApplyOp_+= this->blockSize_;
      OPT::Apply(*this->Op_,*this->V_,*this->KV_);
    }
  }



  template <class ScalarType, class MV, class OP>
  void TraceMin<ScalarType,MV,OP>::harmonicAddToBasis(const Teuchos::RCP<const MV> Delta)
  {
    // V = X - Delta
    MVT::MvAddMv(ONE,*this->X_,-ONE,*Delta,*this->V_);

    // Project out auxVecs
    if(this->numAuxVecs_ > 0)
    {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor lcltimer( *this->timerOrtho_ );
#endif
      this->orthman_->projectMat(*this->V_,this->auxVecs_);
    }

    // Compute KV
    if(this->Op_ != Teuchos::null)
    {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor lcltimer( *this->timerOp_ );
#endif
      this->count_ApplyOp_+= this->blockSize_;

      OPT::Apply(*this->Op_,*this->V_,*this->KV_);
    }

    // Normalize lclKV
    RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > gamma = rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(this->blockSize_,this->blockSize_));
    int rank = this->orthman_->normalizeMat(*this->KV_,gamma);
    // FIXME (mfh 18 Feb 2015) It would make sense to check the rank.
    (void) rank;

    // lclV = lclV/gamma
    Teuchos::SerialDenseSolver<int,ScalarType> SDsolver;
    SDsolver.setMatrix(gamma);
    SDsolver.invert();
    RCP<MV> tempMV = MVT::CloneCopy(*this->V_);
    MVT::MvTimesMatAddMv(ONE,*tempMV,*gamma,ZERO,*this->V_);

    // Compute MV
    if(this->hasM_)
    {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor lcltimer( *this->timerMOp_ );
#endif
      this->count_ApplyM_+= this->blockSize_;

      OPT::Apply(*this->MOp_,*this->V_,*this->MV_);
    }
  }

}} // End of namespace Anasazi

#endif

// End of file AnasaziTraceMin.hpp