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 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_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;
    
    int nvecs = MVT::GetNumberVecs(X);
    
    // Make sure that normvec has enough entries to hold the norms
    // of all the columns of X.  std::vector<T>::size_type is
    // unsigned, so do the appropriate cast to avoid signed/unsigned
    // comparisons that trigger compiler warnings.
    if (normvec.size() < static_cast<size_t>(nvecs))
      normvec.resize (nvecs);
    
    if (!_hasOp) {
      // X == MX, since the operator M is the identity.
      MX = Teuchos::rcp(&X, false);
      MVT::MvNorm(X, normvec);
    }
    else {
      // The caller didn't give us a previously computed MX, so
      // apply the operator.  We assign to MX only after applying
      // the operator, so that if the application fails, MX won't be
      // modified.
      if(MX == Teuchos::null) {
        Teuchos::RCP<MV> tempVec = MVT::Clone(X,nvecs);
        OPT::Apply(*_Op,X,*tempVec);
        _OpCounter += nvecs;
        MX = tempVec;
      }
      else {
        // The caller gave us a previously computed MX.  Make sure
        // that it has at least as many columns as X.
        const int numColsMX = MVT::GetNumberVecs(*MX);
        TEUCHOS_TEST_FOR_EXCEPTION(numColsMX < nvecs, std::invalid_argument,
                           "MatOrthoManager::norm(X, MX, normvec): "
                           "MX has fewer columns than X: "
                           "MX has " << numColsMX << " columns, "
                           "and X has " << nvecs << " columns.");
      }
      
      std::vector<ScalarType> dotvec(nvecs);
      MVT::MvDot(X,*MX,dotvec);
      for (int i=0; i<nvecs; i++) {
        normvec[i] = MT::squareroot( SCT::magnitude(dotvec[i]) );
      }
    }
  }

  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