Ifpack2_RILUK_decl.hpp

/*@HEADER
// ***********************************************************************
// 
//       Ifpack2: Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2002) 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 IFPACK2_CRSRILUK_DECL_HPP
#define IFPACK2_CRSRILUK_DECL_HPP

#include "Ifpack2_ConfigDefs.hpp"
#include "Ifpack2_ScalingType.hpp"
#include "Ifpack2_IlukGraph.hpp"
#include "Ifpack2_Preconditioner.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_MultiVector.hpp"

#include "Teuchos_RefCountPtr.hpp"

namespace Teuchos {
  class ParameterList;
}

namespace Ifpack2 {


template<class MatrixType>
class RILUK: public virtual Ifpack2::Preconditioner<typename MatrixType::scalar_type,typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> {
      
 public:
  typedef typename MatrixType::scalar_type Scalar;
  typedef typename MatrixType::local_ordinal_type LocalOrdinal;
  typedef typename MatrixType::global_ordinal_type GlobalOrdinal;
  typedef typename MatrixType::node_type Node;
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;


  RILUK(const Teuchos::RCP<const MatrixType>& A_in);
  
 private:
  RILUK(const RILUK<MatrixType> & src);

 public:
  virtual ~RILUK();

  void SetRelaxValue( magnitudeType RelaxValue) {RelaxValue_ = RelaxValue;}

  void SetAbsoluteThreshold( magnitudeType Athresh) {Athresh_ = Athresh;}

  void SetRelativeThreshold( magnitudeType Rthresh) {Rthresh_ = Rthresh;}

  void SetOverlapMode( Tpetra::CombineMode OverlapMode) {OverlapMode_ = OverlapMode;}


  void setParameters(const Teuchos::ParameterList& parameterlist);

  void initialize();
  bool isInitialized() const {return isInitialized_;}
  int getNumInitialize() const {return numInitialize_;}


  void compute();

  bool isComputed() const {return(Factored_);}
  
  int getNumCompute() const {return numCompute_;}
  int getNumApply() const {return numApply_;}

  double getInitializeTime() const {return -1;}
  double getComputeTime() const {return -1;}
  double getApplyTime() const {return -1;}

  // Mathematical functions.
  
  

  void apply(
      const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
            Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
            Teuchos::ETransp mode = Teuchos::NO_TRANS,
               Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
               Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const;



  int Multiply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
                     Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
               Teuchos::ETransp mode = Teuchos::NO_TRANS) const;


  magnitudeType computeCondEst(Teuchos::ETransp mode) const;
  magnitudeType computeCondEst(CondestType CT = Ifpack2::Cheap,
                               LocalOrdinal MaxIters = 1550,
                               magnitudeType Tol = 1e-9,
                               const Teuchos::Ptr<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &Matrix = Teuchos::null)
  {
    std::cerr << "Warning, Ifpack2::RILUK::computeCondEst currently does not use MaxIters/Tol/etc arguments..." << std::endl;
    return computeCondEst(Teuchos::NO_TRANS);
  }

  magnitudeType getCondEst() const {return Condest_;}

  Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > getMatrix() const
  {
    return A_;
  }

  // Attribute access functions
  
  magnitudeType GetRelaxValue() const {return RelaxValue_;}

  magnitudeType getAbsoluteThreshold() const {return Athresh_;}

  magnitudeType getRelativeThreshold() const {return Rthresh_;}

  int getLevelOfFill() const { return LevelOfFill_; }

  Tpetra::CombineMode getOverlapMode() {return OverlapMode_;}
 
  int getGlobalNumEntries() const {return(getL().getGlobalNumEntries()+getU().getGlobalNumEntries());}
 
  const Teuchos::RCP<Ifpack2::IlukGraph<LocalOrdinal,GlobalOrdinal,Node> >& getGraph() const {return(Graph_);}

  const MatrixType& getL() const {return(*L_);}
    
  const Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & getD() const {return(*D_);}

  const MatrixType& getU() const {return(*U_);}


  const Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& getDomainMap() const
  { return Graph_->getL_Graph()->getDomainMap(); }

  const Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& getRangeMap() const
  { return Graph_->getU_Graph()->getRangeMap(); }


 protected:
  void setFactored(bool Flag) {Factored_ = Flag;}
  void setInitialized(bool Flag) {isInitialized_ = Flag;}
  bool isAllocated() const {return(isAllocated_);}
  void setAllocated(bool Flag) {isAllocated_ = Flag;}
  
 private:
  
  
  void allocate_L_and_U();
  void initAllValues(const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> & overlapA);
  void generateXY(Teuchos::ETransp mode, 
     const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Xin,
     const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Yin,
     Teuchos::RCP<const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Xout, 
     Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Yout) const;
  bool isOverlapped_;
  Teuchos::RCP<Ifpack2::IlukGraph<LocalOrdinal,GlobalOrdinal,Node> > Graph_;
  const Teuchos::RCP<const MatrixType> A_;
  Teuchos::RCP<MatrixType> L_;
  Teuchos::RCP<MatrixType> U_;
  Teuchos::RCP<Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > D_;
  bool UseTranspose_;

  int LevelOfFill_;
  int LevelOfOverlap_;

  int NumMyDiagonals_;
  bool isAllocated_;
  bool isInitialized_;
  mutable int numInitialize_;
  mutable int numCompute_;
  mutable int numApply_;
  bool Factored_;
  magnitudeType RelaxValue_;
  magnitudeType Athresh_;
  magnitudeType Rthresh_;
  mutable magnitudeType Condest_;

  mutable Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > OverlapX_;
  mutable Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > OverlapY_;
  Tpetra::CombineMode OverlapMode_;
};

}//namespace Ifpack2

#endif /* IFPACK2_CRSRILUK_DECL_HPP */