Ifpack_AdditiveSchwarz.h

#ifndef IFPACK_ADDITIVESCHWARZ_H
#define IFPACK_ADDITIVESCHWARZ_H

#include "Ifpack_ConfigDefs.h"
#if defined(HAVE_IFPACK_TEUCHOS)
#include "Ifpack_Preconditioner.h"
#include "Teuchos_ParameterList.hpp"
#include "Ifpack_ConfigDefs.h"
#include "Ifpack_Preconditioner.h"
#include "Ifpack_Reordering.h"
#include "Ifpack_RCMReordering.h"
#include "Ifpack_METISReordering.h"
#include "Ifpack_LocalFilter.h"
#include "Ifpack_SingletonFilter.h"
#include "Ifpack_ReorderFilter.h"
#include "Ifpack_Utils.h"
#include "Ifpack_OverlappingRowMatrix.h"
#include "Epetra_CombineMode.h"
#include "Epetra_MultiVector.h"
#include "Epetra_Map.h"
#include "Epetra_Comm.h"
#include "Epetra_Time.h"
#include "Epetra_LinearProblem.h"
#include "Epetra_RowMatrix.h"
#include "Epetra_CrsMatrix.h"
#include "Teuchos_ParameterList.hpp"


template<typename T>
class Ifpack_AdditiveSchwarz : public virtual Ifpack_Preconditioner {
      
public:



  Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix,
             int OverlapLevel = 0);
  
  virtual ~Ifpack_AdditiveSchwarz();



    virtual int SetUseTranspose(bool UseTranspose);
  


    virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;


    virtual int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;

    virtual double NormInf() const;
  

    virtual const char * Label() const;

    virtual bool UseTranspose() const;

    virtual bool HasNormInf() const;

    virtual const Epetra_Comm & Comm() const;

    virtual const Epetra_Map & OperatorDomainMap() const;

    virtual const Epetra_Map & OperatorRangeMap() const;

  virtual bool IsInitialized() const
  {
    return(IsInitialized_);
  }

  virtual bool IsComputed() const
  {
    return(IsComputed_);
  }


  virtual int SetParameters(Teuchos::ParameterList& List);

  // @}

  // @{ Query methods
  
  virtual int Initialize();

  virtual int Compute();

  virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
                         const int MaxIters = 1550,
                         const double Tol = 1e-9,
             Epetra_RowMatrix* Matrix = 0);

  virtual double Condest() const
  {
    return(Condest_);
  }

  virtual const Epetra_RowMatrix& Matrix() const
  {
    return(*Matrix_);
  }

  virtual bool IsOverlapping() const
  {
    return(IsOverlapping_);
  }

  virtual std::ostream& Print(std::ostream&) const;
  
  virtual const T* Inverse() const
  {
    return(Inverse_);
  }

  virtual int NumInitialize() const
  {
    return(NumInitialize_);
  }

  virtual int NumCompute() const
  {
    return(NumCompute_);
  }

  virtual int NumApplyInverse() const
  {
    return(NumApplyInverse_);
  }

  virtual double InitializeTime() const
  {
    return(InitializeTime_);
  }

  virtual double ComputeTime() const
  {
    return(ComputeTime_);
  }

  virtual double ApplyInverseTime() const
  {
    return(ApplyInverseTime_);
  }

  virtual double InitializeFlops() const
  {
    return(InitializeFlops_);
  }

  virtual double ComputeFlops() const
  {
    return(ComputeFlops_);
  }

  virtual double ApplyInverseFlops() const
  {
    return(ApplyInverseFlops_);
  }

  virtual int OverlapLevel() const
  {
    return(OverlapLevel_);
  }

  virtual const Teuchos::ParameterList& List() const
  {
    return(List_);
  }

protected:

  // @}

  // @{ Internal merhods.
  
  Ifpack_AdditiveSchwarz(const Ifpack_AdditiveSchwarz& RHS)
  { }

  int Setup();
  
  void Destroy();

  // @}

  // @{ Internal data.
  
  const Epetra_RowMatrix* Matrix_;
  Ifpack_OverlappingRowMatrix* OverlappingMatrix_;
  Ifpack_LocalFilter* LocalizedMatrix_;
  string Label_;
  bool IsInitialized_;
  bool IsComputed_;
  T* Inverse_;
  bool UseTranspose_;
  bool IsOverlapping_;
  int OverlapLevel_;
  Teuchos::ParameterList List_;
  Epetra_CombineMode CombineMode_;
  double Condest_;
  bool ComputeCondest_;
  bool UseReordering_;
  string ReorderingType_;
  Ifpack_Reordering* Reordering_;
  Ifpack_ReorderFilter* ReorderedLocalizedMatrix_;
  bool FilterSingletons_;
  Ifpack_SingletonFilter* SingletonFilter_;
  int NumInitialize_;
  int NumCompute_;
  mutable int NumApplyInverse_;
  double InitializeTime_;
  double ComputeTime_;
  mutable double ApplyInverseTime_;
  double InitializeFlops_;
  double ComputeFlops_;
  mutable double ApplyInverseFlops_;
  Epetra_Time* Time_;

}; // class Ifpack_AdditiveSchwarz<T>

//==============================================================================
template<typename T>
Ifpack_AdditiveSchwarz<T>::
Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix,
               int OverlapLevel) :
  Matrix_(Matrix),
  OverlappingMatrix_(0),
  LocalizedMatrix_(0),
  IsInitialized_(false),
  IsComputed_(false),
  Inverse_(0),
  UseTranspose_(false),
  IsOverlapping_(false),
  OverlapLevel_(OverlapLevel),
  CombineMode_(Zero),
  Condest_(-1.0),
  ComputeCondest_(true),
  UseReordering_(false),
  ReorderingType_("none"),
  Reordering_(0),
  ReorderedLocalizedMatrix_(0),
  FilterSingletons_(false),
  SingletonFilter_(0),
  NumInitialize_(0),
  NumCompute_(0),
  NumApplyInverse_(0),
  InitializeTime_(0.0),
  ComputeTime_(0.0),
  ApplyInverseTime_(0.0),
  InitializeFlops_(0.0),
  ComputeFlops_(0.0),
  ApplyInverseFlops_(0.0),
  Time_(0)
{
  if (Matrix_->Comm().NumProc() == 1)
    OverlapLevel_ = 0;

  if ((OverlapLevel_ != 0) && (Matrix_->Comm().NumProc() > 1))
    IsOverlapping_ = true;
  // Sets parameters to default values
  Teuchos::ParameterList List;
  SetParameters(List);
}

//==============================================================================
template<typename T>
Ifpack_AdditiveSchwarz<T>::~Ifpack_AdditiveSchwarz()
{
  Destroy();
}

//==============================================================================
template<typename T>
void Ifpack_AdditiveSchwarz<T>::Destroy() 
{
  if (OverlappingMatrix_)
    delete OverlappingMatrix_;
  OverlappingMatrix_ = 0;

  if (Inverse_)
    delete Inverse_;
  Inverse_ = 0;

  if (LocalizedMatrix_)
    delete LocalizedMatrix_;
  LocalizedMatrix_ = 0;

  if (ReorderedLocalizedMatrix_)
    delete ReorderedLocalizedMatrix_;
  ReorderedLocalizedMatrix_ = 0;

  if (SingletonFilter_)
    delete SingletonFilter_;
  SingletonFilter_ = 0;

  if (Reordering_)
    delete Reordering_;
  Reordering_ = 0;

  if (Time_)
    delete Time_;
  Time_ = 0;
}

//==============================================================================
template<typename T>
int Ifpack_AdditiveSchwarz<T>::Setup()
{

  Epetra_RowMatrix* MatrixPtr;

  if (OverlappingMatrix_)
    LocalizedMatrix_ = new Ifpack_LocalFilter(OverlappingMatrix_);
  else
    LocalizedMatrix_ = new Ifpack_LocalFilter(Matrix_);

  if (LocalizedMatrix_ == 0)
    IFPACK_CHK_ERR(-5);

  // users may want to skip singleton check
  if (FilterSingletons_) {
    SingletonFilter_ = new Ifpack_SingletonFilter(LocalizedMatrix_);
    MatrixPtr = SingletonFilter_;
  }
  else
    MatrixPtr = LocalizedMatrix_;

  if (UseReordering_) {

    // create reordeing and compute it
    if (ReorderingType_ == "rcm")
      Reordering_ = new Ifpack_RCMReordering();
    else if (ReorderingType_ == "metis")
      Reordering_ = new Ifpack_METISReordering();
    else {
      cerr << "reordering type not correct (" << ReorderingType_ << ")" << endl;
      exit(EXIT_FAILURE);
    }
    if (Reordering_ == 0) IFPACK_CHK_ERR(-5);

    IFPACK_CHK_ERR(Reordering_->SetParameters(List_));
    IFPACK_CHK_ERR(Reordering_->Compute(*MatrixPtr));

    // now create reordered localized matrix
    ReorderedLocalizedMatrix_ = 
      new Ifpack_ReorderFilter(MatrixPtr,Reordering_);

    if (ReorderedLocalizedMatrix_ == 0) IFPACK_CHK_ERR(-5);

    MatrixPtr = ReorderedLocalizedMatrix_;
  }

  Inverse_ = new T(MatrixPtr);

  if (Inverse_ == 0)
    IFPACK_CHK_ERR(-5);

  return(0);
}

//==============================================================================
template<typename T>
int Ifpack_AdditiveSchwarz<T>::SetParameters(Teuchos::ParameterList& List)
{
 
  // compute the condition number each time Compute() is invoked.
  ComputeCondest_ = List.get("schwarz: compute condest", ComputeCondest_);
  // combine mode
  if( Teuchos::ParameterEntry *combineModeEntry = List.getEntryPtr("schwarz: combine mode") )
  {
    if( typeid(std::string) == combineModeEntry->getAny().type() )
    {
      std::string mode = List.get("schwarz: combine mode", "Add");
      if (mode == "Add")
        CombineMode_ = Add;
      else if (mode == "Zero")
        CombineMode_ = Zero;
      else if (mode == "Insert")
        CombineMode_ = Insert;
      else if (mode == "InsertAdd")
        CombineMode_ = InsertAdd;
      else if (mode == "Average")
        CombineMode_ = Average;
      else if (mode == "AbsMax")
        CombineMode_ = AbsMax;
      else
      {
        TEST_FOR_EXCEPTION(
          true,std::logic_error
          ,"Error, The (Epetra) combine mode of \""<<mode<<"\" is not valid!  Only the values"
          " \"Add\", \"Zero\", \"Insert\", \"InsertAdd\", \"Average\", and \"AbsMax\" are accepted!"
          );
      }
    }
    else if ( typeid(Epetra_CombineMode) == combineModeEntry->getAny().type() )
    {
      CombineMode_ = Teuchos::any_cast<Epetra_CombineMode>(combineModeEntry->getAny());
    }
    else
    {
      // Throw exception with good error message!
      Teuchos::getParameter<std::string>(List,"schwarz: combine mode");
    }
  }
  else
  {
    // Make the default be a string to be consistent with the valid parameters!
    List.get("schwarz: combine mode","Zero");
  }
  // type of reordering
  ReorderingType_ = List.get("schwarz: reordering type", ReorderingType_);
  if (ReorderingType_ == "none")
    UseReordering_ = false;
  else 
    UseReordering_ = true;
  // if true, filter singletons. NOTE: the filtered matrix can still have
  // singletons! A simple example: upper triangular matrix, if I remove
  // the lower node, I still get a matrix with a singleton! However, filter
  // singletons should help for PDE problems with Dirichlet BCs.
  FilterSingletons_ = List.get("schwarz: filter singletons", FilterSingletons_);

  // This copy may be needed by Amesos or other preconditioners.
  List_ = List;

  return(0);
}

//==============================================================================
template<typename T>
int Ifpack_AdditiveSchwarz<T>::Initialize()
{
  IsInitialized_ = false;
  IsComputed_ = false; // values required
  Condest_ = -1.0; // zero-out condest

  Destroy();

  if (Time_ == 0)
    Time_ = new Epetra_Time(Comm());

  Time_->ResetStartTime();

  // compute the overlapping matrix if necessary
  if (IsOverlapping_) {
    OverlappingMatrix_ = 
      new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_);
    if (OverlappingMatrix_ == 0)
      IFPACK_CHK_ERR(-5);
  }

  IFPACK_CHK_ERR(Setup());

  if (Inverse_ == 0)
    IFPACK_CHK_ERR(-5);

  if (LocalizedMatrix_ == 0)
    IFPACK_CHK_ERR(-5);

  IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose()));
  IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
  IFPACK_CHK_ERR(Inverse_->Initialize());

  // Label is for Aztec-like solvers
  Label_ = "Ifpack_AdditiveSchwarz, ";
  if (UseTranspose())
    Label_ += ", transp";
  Label_ += ", ov = " + Ifpack_toString(OverlapLevel_)
    + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'";

  IsInitialized_ = true;
  ++NumInitialize_;
  InitializeTime_ += Time_->ElapsedTime();

  // count flops by summing up all the InitializeFlops() in each
  // Inverse. Each Inverse() can only give its flops -- it acts on one
  // process only
  double partial = Inverse_->InitializeFlops();
  double total;
  Comm().SumAll(&partial, &total, 1);
  InitializeFlops_ += total;

  return(0);
}

//==============================================================================
template<typename T>
int Ifpack_AdditiveSchwarz<T>::Compute()
{

  if (IsInitialized() == false)
    IFPACK_CHK_ERR(Initialize());

  Time_->ResetStartTime();
  IsComputed_ = false;
  Condest_ = -1.0;
  
  IFPACK_CHK_ERR(Inverse_->Compute());

  IsComputed_ = true; // need this here for Condest(Ifpack_Cheap)
  ++NumCompute_;
  ComputeTime_ += Time_->ElapsedTime();

  // sum up flops
  double partial = Inverse_->ComputeFlops();
  double total;
  Comm().SumAll(&partial, &total, 1);
  ComputeFlops_ += total;

  // reset the Label
  string R = "";
  if (UseReordering_)
    R = ReorderingType_ + " reord, ";

  if (ComputeCondest_)
    Condest(Ifpack_Cheap);
  
  // add Condest() to label
  Label_ = "Ifpack_AdditiveSchwarz, ov = " + Ifpack_toString(OverlapLevel_)
    + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'"
    + "\n\t\t***** " + R + "Condition number estimate = "
    + Ifpack_toString(Condest());

  return(0);
}

//==============================================================================
template<typename T>
int Ifpack_AdditiveSchwarz<T>::SetUseTranspose(bool UseTranspose)
{
  // store the flag -- it will be set in Initialize() if Inverse_ does not
  // exist.
  UseTranspose_ = UseTranspose;

  // If Inverse_ exists, pass it right now.
  if (Inverse_)
    IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose));
  return(0);
}

//==============================================================================
template<typename T>
int Ifpack_AdditiveSchwarz<T>::
Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{
  IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
  return(0);
}

//==============================================================================
template<typename T>
double Ifpack_AdditiveSchwarz<T>::NormInf() const
{
  return(-1.0);
}

//==============================================================================
template<typename T>
const char * Ifpack_AdditiveSchwarz<T>::Label() const
{
  return(Label_.c_str());
}

//==============================================================================
template<typename T>
bool Ifpack_AdditiveSchwarz<T>::UseTranspose() const
{
  return(UseTranspose_);
}

//==============================================================================
template<typename T>
bool Ifpack_AdditiveSchwarz<T>::HasNormInf() const
{
  return(false);
}

//==============================================================================
template<typename T>
const Epetra_Comm & Ifpack_AdditiveSchwarz<T>::Comm() const
{
  return(Matrix_->Comm());
}

//==============================================================================
template<typename T>
const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorDomainMap() const
{
  return(Matrix_->OperatorDomainMap());
}

//==============================================================================
template<typename T>
const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorRangeMap() const
{
  return(Matrix_->OperatorRangeMap());
}

//==============================================================================
template<typename T>
int Ifpack_AdditiveSchwarz<T>::
ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{
  // compute the preconditioner is not done by the user
  if (!IsComputed())
    IFPACK_CHK_ERR(-3);
  
  int NumVectors = X.NumVectors();

  if (NumVectors != Y.NumVectors())
    IFPACK_CHK_ERR(-2); // wrong input

  Time_->ResetStartTime();

  Epetra_MultiVector* OverlappingX;
  Epetra_MultiVector* OverlappingY;
  Epetra_MultiVector* Xtmp = 0;

  // for flop count, see bottom of this function
  double pre_partial = Inverse_->ApplyInverseFlops();
  double pre_total;
  Comm().SumAll(&pre_partial, &pre_total, 1);

  // process overlap, may need to create vectors and import data
  if (IsOverlapping()) {
    OverlappingX = new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
                      X.NumVectors());
    OverlappingY = new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
                                          Y.NumVectors());
    if (OverlappingY == 0) IFPACK_CHK_ERR(-5);

    OverlappingY->PutScalar(0.0);
    OverlappingX->PutScalar(0.0);
    IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(X,*OverlappingX,Insert));
    // FIXME: this will not work with overlapping and non-zero starting
    // solutions. The same for other cases below.
    // IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(Y,*OverlappingY,Insert));
  }
  else {
    Xtmp = new Epetra_MultiVector(X);
    OverlappingX = (Epetra_MultiVector*)Xtmp;
    OverlappingY = &Y;
  }

  if (FilterSingletons_) {
    // process singleton filter
    Epetra_MultiVector ReducedX(SingletonFilter_->Map(),NumVectors);
    Epetra_MultiVector ReducedY(SingletonFilter_->Map(),NumVectors);
    IFPACK_CHK_ERR(SingletonFilter_->SolveSingletons(*OverlappingX,*OverlappingY));
    IFPACK_CHK_ERR(SingletonFilter_->CreateReducedRHS(*OverlappingY,*OverlappingX,ReducedX));

    // process reordering
    if (!UseReordering_) {
      IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReducedX,ReducedY));
    }
    else {
      Epetra_MultiVector ReorderedX(ReducedX);
      Epetra_MultiVector ReorderedY(ReducedY);
      IFPACK_CHK_ERR(Reordering_->P(ReducedX,ReorderedX));
      IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
      IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,ReducedY));
    }

    // finish up with singletons
    IFPACK_CHK_ERR(SingletonFilter_->UpdateLHS(ReducedY,*OverlappingY));
  }
  else {
    // process reordering
    if (!UseReordering_) {
      IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX,*OverlappingY));
    }
    else {
      Epetra_MultiVector ReorderedX(*OverlappingX);
      Epetra_MultiVector ReorderedY(*OverlappingY);
      IFPACK_CHK_ERR(Reordering_->P(*OverlappingX,ReorderedX));
      IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
      IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,*OverlappingY));
    }
  }

  if (IsOverlapping()) {
    IFPACK_CHK_ERR(OverlappingMatrix_->ExportMultiVector(*OverlappingY,Y,
                                                        CombineMode_));

    delete OverlappingX;
    delete OverlappingY;
  }

  if (Xtmp)
    delete Xtmp;

  // add flops. Note the we only have to add the newly counted
  // flops -- and each Inverse returns the cumulative sum
  double partial = Inverse_->ApplyInverseFlops();
  double total;
  Comm().SumAll(&partial, &total, 1);
  ApplyInverseFlops_ += total - pre_total;

  // FIXME: right now I am skipping the overlap and singletons
  ++NumApplyInverse_;
  ApplyInverseTime_ += Time_->ElapsedTime();

  return(0);
 
}

//==============================================================================
template<typename T>
std::ostream& Ifpack_AdditiveSchwarz<T>::
Print(std::ostream& os) const
{
  double IF = InitializeFlops();
  double CF = ComputeFlops();
  double AF = ApplyInverseFlops();

  double IFT = 0.0, CFT = 0.0, AFT = 0.0;
  if (InitializeTime() != 0.0)
    IFT = IF / InitializeTime();
  if (ComputeTime() != 0.0)
    CFT = CF / ComputeTime();
  if (ApplyInverseTime() != 0.0)
    AFT = AF / ApplyInverseTime();

  if (Matrix().Comm().MyPID())
    return(os);

  os << endl;
  os << "================================================================================" << endl;
  os << "Ifpack_AdditiveSchwarz, overlap level = " << OverlapLevel_ << endl;
  if (CombineMode_ == Insert)
    os << "Combine mode                          = Insert" << endl;
  else if (CombineMode_ == Add)
    os << "Combine mode                          = Add" << endl;
  else if (CombineMode_ == Zero)
    os << "Combine mode                          = Zero" << endl;
  else if (CombineMode_ == Average)
    os << "Combine mode                          = Average" << endl;
  else if (CombineMode_ == AbsMax)
    os << "Combine mode                          = AbsMax" << endl;

  os << "Condition number estimate             = " << Condest_ << endl;
  os << "Global number of rows                 = " << Matrix_->NumGlobalRows() << endl;

  os << endl;
  os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
  os << "-----           -------   --------------       ------------     --------" << endl;
  os << "Initialize()    "   << std::setw(5) << NumInitialize()
     << "  " << std::setw(15) << InitializeTime() 
     << "  " << std::setw(15) << 1.0e-6 * IF 
     << "  " << std::setw(15) << 1.0e-6 * IFT << endl;
  os << "Compute()       "   << std::setw(5) << NumCompute() 
     << "  " << std::setw(15) << ComputeTime()
     << "  " << std::setw(15) << 1.0e-6 * CF
     << "  " << std::setw(15) << 1.0e-6 * CFT << endl;
  os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
     << "  " << std::setw(15) << ApplyInverseTime()
     << "  " << std::setw(15) << 1.0e-6 * AF
     << "  " << std::setw(15) << 1.0e-6 * AFT << endl;
  os << "================================================================================" << endl;
  os << endl;

  return(os);
}

#include "Ifpack_Condest.h"
//==============================================================================
template<typename T>
double Ifpack_AdditiveSchwarz<T>::
Condest(const Ifpack_CondestType CT, const int MaxIters, 
        const double Tol, Epetra_RowMatrix* Matrix)
{
  if (!IsComputed()) // cannot compute right now
    return(-1.0);

  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);

  return(Condest_);
}

#endif // HAVE_IFPACK_TEUCHOS
#endif // IFPACK_ADDITIVESCHWARZ_H

---