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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// 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> {
  protected:
    Teuchos::RefCountPtr<const OP> _Op;
    bool _hasOp;

  public:


    MatOrthoManager(Teuchos::RefCountPtr<const OP> Op = Teuchos::null) : _Op(Op), _hasOp(Op!=Teuchos::null) {};

    virtual ~MatOrthoManager() {};



    void setOp( Teuchos::RefCountPtr<const OP> Op ) { 
      _Op = Op; 
      _hasOp = (_Op != Teuchos::null);
    };

    Teuchos::RefCountPtr<const OP> getOp() const { return _Op; } 





    void 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::RefCountPtr<const MV> P,Q;
      Teuchos::RefCountPtr<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);
          P = R;
          Q = Teuchos::rcp( &Y, false );
        }
        else {
          P = Teuchos::rcp( &X, false );
          R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
          OPT::Apply(*_Op,Y,*R);
          Q = R;
        }
      }
      else {
        P = Teuchos::rcp( &X, false );
        Q = Teuchos::rcp( &Y, false );
      }
      
      MVT::MvTransMv(SCT::one(),*P,*Q,Z);
    }

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

      Teuchos::RefCountPtr<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);
      }
    }

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

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

      typedef Teuchos::ScalarTraits<ScalarType> SCT;
      typedef MultiVecTraits<ScalarType,MV>     MVT;
      typedef OperatorTraits<ScalarType,MV,OP>  OPT;
      
      if (!_hasOp) {
        MX = Teuchos::rcp(&X,false);
      }
      else if (MX == Teuchos::null) {
        Teuchos::RefCountPtr<MV> R = MVT::Clone(X,MVT::GetNumberVecs(X));
        OPT::Apply(*_Op,X,*R);
        MX = R;
      }

      Teuchos::SerialDenseMatrix<int,ScalarType> z(1,1);
      Teuchos::RefCountPtr<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] = SCT::magnitude( SCT::squareroot( z(0,0) ) );
      }
    }


    virtual void project ( MV &X, Teuchos::RefCountPtr<MV> MX, 
                           Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
                           Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q) const = 0;


    
    virtual void project ( MV &X, 
                           Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
                           Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q) const {
      project(X,Teuchos::null,C,Q);
    }

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


    virtual int normalize ( MV &X, Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
      return normalize(X,Teuchos::null,B);
    }


    virtual int projectAndNormalize ( MV &X, Teuchos::RefCountPtr<MV> MX, 
                                      Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
                                      Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
                                      Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q ) const = 0;

    virtual int projectAndNormalize ( MV &X, 
                                      Teuchos::Array<Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
                                      Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
                                      Teuchos::Array<Teuchos::RefCountPtr<const MV> > Q ) const {
      return projectAndNormalize(X,Teuchos::null,C,B,Q);
    }




    virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
    orthonormError(const MV &X) const {
      return orthonormError(X,Teuchos::null);
    }

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

    virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
    orthogError(const MV &X1, const MV &X2) const {
      return orthogError(X1,Teuchos::null,X2);
    }

    virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
    orthogError(const MV &X1, Teuchos::RefCountPtr<const MV> MX1, const MV &X2) const = 0;


  };

} // end of Anasazi namespace


#endif

// end of file AnasaziMatOrthoManager.hpp

---