AnasaziMatOrthoManager.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_MATORTHOMANAGER_HPP
#define ANASAZI_MATORTHOMANAGER_HPP

#include "AnasaziConfigDefs.hpp"
#include "AnasaziTypes.hpp"
#include "AnasaziOrthoManager.hpp"
#include "AnasaziMultiVecTraits.hpp"
#include "AnasaziOperatorTraits.hpp"

namespace Anasazi {

  template <class ScalarType, class MV, class OP>
  class MatOrthoManager : public OrthoManager<ScalarType,MV> {
  public:


    MatOrthoManager(Teuchos::RCP<const OP> Op = Teuchos::null);

    virtual ~MatOrthoManager() {};



    virtual void setOp( Teuchos::RCP<const OP> Op );

    virtual Teuchos::RCP<const OP> getOp() const;


    int getOpCounter() const;


    void resetOpCounter();




    void innerProdMat(
          const MV& X, const MV& Y,
          Teuchos::SerialDenseMatrix<int,ScalarType>& Z,
          Teuchos::RCP<const MV> MX = Teuchos::null,
          Teuchos::RCP<const MV> MY = Teuchos::null
        ) const;

    void normMat(
          const MV& X,
          std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
          Teuchos::RCP<const MV> MX = Teuchos::null
        ) const;

    virtual void projectMat (
          MV &X,
          Teuchos::Array<Teuchos::RCP<const MV> >  Q,
          Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
              = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
          Teuchos::RCP<MV> MX                                                          = Teuchos::null,
          Teuchos::Array<Teuchos::RCP<const MV> > MQ
              = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
        ) const = 0;

    virtual int normalizeMat (
          MV &X,
          Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
          Teuchos::RCP<MV> MX                                         = Teuchos::null
        ) const = 0;


    virtual int projectAndNormalizeMat (
          MV &X,
          Teuchos::Array<Teuchos::RCP<const MV> >  Q,
          Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
              = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
          Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B                  = Teuchos::null,
          Teuchos::RCP<MV> MX                                                          = Teuchos::null,
          Teuchos::Array<Teuchos::RCP<const MV> > MQ
              = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
        ) const = 0;

    virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
    orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const = 0;

    virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
    orthogErrorMat(
          const MV &X,
          const MV &Y,
          Teuchos::RCP<const MV> MX = Teuchos::null,
          Teuchos::RCP<const MV> MY = Teuchos::null
        ) const = 0;




    void innerProd( const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const;

    void norm( const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec ) const;

    void project (
          MV &X,
          Teuchos::Array<Teuchos::RCP<const MV> > Q,
          Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
              = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null))
        ) const;

    int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null) const;

    int projectAndNormalize (
          MV &X,
          Teuchos::Array<Teuchos::RCP<const MV> > Q,
          Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
              = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
          Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null
        ) const;

    typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
    orthonormError(const MV &X) const;

    typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
    orthogError(const MV &X1, const MV &X2) const;


  protected:
    Teuchos::RCP<const OP> _Op;
    bool _hasOp;
    mutable int _OpCounter;

  };

  template <class ScalarType, class MV, class OP>
  MatOrthoManager<ScalarType,MV,OP>::MatOrthoManager(Teuchos::RCP<const OP> Op)
      : _Op(Op), _hasOp(Op!=Teuchos::null), _OpCounter(0) {}

  template <class ScalarType, class MV, class OP>
  void MatOrthoManager<ScalarType,MV,OP>::setOp( Teuchos::RCP<const OP> Op )
  {
    _Op = Op;
    _hasOp = (_Op != Teuchos::null);
  }

  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<const OP> MatOrthoManager<ScalarType,MV,OP>::getOp() const
  {
    return _Op;
  }

  template <class ScalarType, class MV, class OP>
  int MatOrthoManager<ScalarType,MV,OP>::getOpCounter() const
  {
    return _OpCounter;
  }

  template <class ScalarType, class MV, class OP>
  void MatOrthoManager<ScalarType,MV,OP>::resetOpCounter()
  {
    _OpCounter = 0;
  }

  template <class ScalarType, class MV, class OP>
  void MatOrthoManager<ScalarType,MV,OP>::innerProd(
      const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const
  {
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef MultiVecTraits<ScalarType,MV>     MVT;
    typedef OperatorTraits<ScalarType,MV,OP>  OPT;

    Teuchos::RCP<const MV> P,Q;
    Teuchos::RCP<MV> R;

    if (_hasOp) {
      // attempt to minimize the amount of work in applying
      if ( MVT::GetNumberVecs(X) < MVT::GetNumberVecs(Y) ) {
        R = MVT::Clone(X,MVT::GetNumberVecs(X));
        OPT::Apply(*_Op,X,*R);
        _OpCounter += MVT::GetNumberVecs(X);
        P = R;
        Q = Teuchos::rcpFromRef(Y);
      }
      else {
        P = Teuchos::rcpFromRef(X);
        R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
        OPT::Apply(*_Op,Y,*R);
        _OpCounter += MVT::GetNumberVecs(Y);
        Q = R;
      }
    }
    else {
      P = Teuchos::rcpFromRef(X);
      Q = Teuchos::rcpFromRef(Y);
    }

    MVT::MvTransMv(SCT::one(),*P,*Q,Z);
  }

  template <class ScalarType, class MV, class OP>
  void MatOrthoManager<ScalarType,MV,OP>::innerProdMat(
      const MV& X, const MV& Y, Teuchos::SerialDenseMatrix<int,ScalarType>& Z, Teuchos::RCP<const MV> MX, Teuchos::RCP<const MV> MY) const
  {
    (void) MX;
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef MultiVecTraits<ScalarType,MV>     MVT;
    // typedef OperatorTraits<ScalarType,MV,OP>  OPT; // unused

    Teuchos::RCP<MV> P,Q;

    if ( MY == Teuchos::null ) {
      innerProd(X,Y,Z);
    }
    else if ( _hasOp ) {
      // the user has done the matrix vector for us
      MVT::MvTransMv(SCT::one(),X,*MY,Z);
    }
    else {
      // there is no matrix vector
      MVT::MvTransMv(SCT::one(),X,Y,Z);
    }
#ifdef TEUCHOS_DEBUG
    for (int j=0; j<Z.numCols(); j++) {
      for (int i=0; i<Z.numRows(); i++) {
        TEUCHOS_TEST_FOR_EXCEPTION(SCT::isnaninf(Z(i,j)), std::logic_error,
            "Anasazi::MatOrthoManager::innerProdMat(): detected NaN/inf.");
      }
    }
#endif
  }

  template <class ScalarType, class MV, class OP>
  void MatOrthoManager<ScalarType,MV,OP>::norm(
      const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec ) const
  {
    this->normMat(X,normvec);
  }

  template <class ScalarType, class MV, class OP>
  void MatOrthoManager<ScalarType,MV,OP>::normMat(
      const MV& X,
      std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &normvec,
      Teuchos::RCP<const MV> MX) const
  {
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef Teuchos::ScalarTraits<typename SCT::magnitudeType> MT;
    typedef MultiVecTraits<ScalarType,MV>     MVT;
    typedef OperatorTraits<ScalarType,MV,OP>  OPT;

    if (!_hasOp) {
      MX = Teuchos::rcpFromRef(X);
    }
    else if (MX == Teuchos::null) {
      Teuchos::RCP<MV> R = MVT::Clone(X,MVT::GetNumberVecs(X));
      OPT::Apply(*_Op,X,*R);
      _OpCounter += MVT::GetNumberVecs(X);
      MX = R;
    }

    Teuchos::SerialDenseMatrix<int,ScalarType> z(1,1);
    Teuchos::RCP<const MV> Xi, MXi;
    std::vector<int> ind(1);
    for (int i=0; i<MVT::GetNumberVecs(X); i++) {
      ind[0] = i;
      Xi = MVT::CloneView(X,ind);
      MXi = MVT::CloneView(*MX,ind);
      MVT::MvTransMv(SCT::one(),*Xi,*MXi,z);
      normvec[i] = MT::squareroot( SCT::magnitude(z(0,0)) );
    }
  }

  template <class ScalarType, class MV, class OP>
  void MatOrthoManager<ScalarType,MV,OP>::project (
        MV &X,
        Teuchos::Array<Teuchos::RCP<const MV> > Q,
        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
      ) const
  {
    this->projectMat(X,Q,C);
  }

  template <class ScalarType, class MV, class OP>
  int MatOrthoManager<ScalarType,MV,OP>::normalize (
      MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const
  {
    return this->normalizeMat(X,B);
  }

  template <class ScalarType, class MV, class OP>
  int MatOrthoManager<ScalarType,MV,OP>::projectAndNormalize (
        MV &X,
        Teuchos::Array<Teuchos::RCP<const MV> > Q,
        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B
      ) const
  {
    return this->projectAndNormalizeMat(X,Q,C,B);
  }

  template <class ScalarType, class MV, class OP>
  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
  MatOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X) const
  {
    return this->orthonormErrorMat(X,Teuchos::null);
  }

  template <class ScalarType, class MV, class OP>
  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
  MatOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, const MV &X2) const
  {
    return this->orthogErrorMat(X1,X2);
  }

} // end of Anasazi namespace


#endif

// end of file AnasaziMatOrthoManager.hpp