Ifpack_AdditiveSchwarz.h

/*
//@HEADER
// ***********************************************************************
//
//       Ifpack: 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.
//
// 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 IFPACK_ADDITIVESCHWARZ_H
#define IFPACK_ADDITIVESCHWARZ_H

#include "Ifpack_ConfigDefs.h"
#include "Ifpack_Preconditioner.h"
#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_NodeFilter.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"
#include "Teuchos_RefCountPtr.hpp"

#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
#include "Ifpack_SubdomainFilter.h"
#include "EpetraExt_Reindex_CrsMatrix.h"
#include "EpetraExt_Reindex_MultiVector.h"
#endif
#ifdef IFPACK_NODE_AWARE_CODE
#include "EpetraExt_OperatorOut.h"
#include "EpetraExt_RowMatrixOut.h"
#include "EpetraExt_BlockMapOut.h"
#endif

#ifdef HAVE_IFPACK_AMESOS
  #include "Ifpack_AMDReordering.h"
#endif



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



  Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix_in,
             int OverlapLevel_in = 0);
  
  virtual ~Ifpack_AdditiveSchwarz() {};



    virtual int SetUseTranspose(bool UseTranspose_in);
  


    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_in = 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();
  
  // @}

  // @{ Internal data.
  
  Teuchos::RefCountPtr<const Epetra_RowMatrix> Matrix_;
  Teuchos::RefCountPtr<Ifpack_OverlappingRowMatrix> OverlappingMatrix_;
/*
  //TODO if we choose to expose the node aware code, i.e., no ifdefs,
  //TODO then we should switch to this definition.
  Teuchos::RefCountPtr<Epetra_RowMatrix> LocalizedMatrix_;
*/
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS

  Teuchos::RCP<Epetra_RowMatrix> LocalizedMatrix_;
  Teuchos::RCP<Epetra_CrsMatrix> SubdomainCrsMatrix_;
  Teuchos::RCP<Epetra_CrsMatrix> ReindexedCrsMatrix_;

  // The following data members are needed when doing ApplyInverse
  // with an Ifpack_SubdomainFilter local matrix
  Teuchos::RCP<Epetra_Map> tempMap_;
  Teuchos::RCP<Epetra_Map> tempDomainMap_;
  Teuchos::RCP<Epetra_Map> tempRangeMap_;
  Teuchos::RCP<EpetraExt::CrsMatrix_Reindex> SubdomainMatrixReindexer_;
  Teuchos::RCP<EpetraExt::MultiVector_Reindex> DomainVectorReindexer_;
  Teuchos::RCP<EpetraExt::MultiVector_Reindex> RangeVectorReindexer_;
  mutable Teuchos::RCP<Epetra_MultiVector> tempX_;
  mutable Teuchos::RCP<Epetra_MultiVector> tempY_;
#else
# ifdef IFPACK_NODE_AWARE_CODE
  Teuchos::RefCountPtr<Ifpack_NodeFilter> LocalizedMatrix_;
# else
  Teuchos::RefCountPtr<Ifpack_LocalFilter> LocalizedMatrix_;
# endif
#endif
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
  // Some data members used for determining the subdomain id (color)
  int MpiRank_;
  int NumMpiProcs_;
  int NumMpiProcsPerSubdomain_;
  int NumSubdomains_;
  int SubdomainId_;
#endif

  string Label_;
  bool IsInitialized_;
  bool IsComputed_;
  bool UseTranspose_;
  bool IsOverlapping_;
  int OverlapLevel_;
  Teuchos::ParameterList List_;
  Epetra_CombineMode CombineMode_;
  double Condest_;
  bool ComputeCondest_;
  bool UseReordering_;
  string ReorderingType_;
  Teuchos::RefCountPtr<Ifpack_Reordering> Reordering_;
  Teuchos::RefCountPtr<Ifpack_ReorderFilter> ReorderedLocalizedMatrix_;
  bool FilterSingletons_;
  Teuchos::RefCountPtr<Ifpack_SingletonFilter> SingletonFilter_;
  int NumInitialize_;
  int NumCompute_;
  mutable int NumApplyInverse_;
  double InitializeTime_;
  double ComputeTime_;
  mutable double ApplyInverseTime_;
  double InitializeFlops_;
  double ComputeFlops_;
  mutable double ApplyInverseFlops_;
  Teuchos::RefCountPtr<Epetra_Time> Time_;
  Teuchos::RefCountPtr<T> Inverse_;
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
  mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
  mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
#endif
# ifdef IFPACK_NODE_AWARE_CODE
  mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
  mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
#endif

}; // class Ifpack_AdditiveSchwarz<T>

//==============================================================================
template<typename T>
Ifpack_AdditiveSchwarz<T>::
Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix_in,
                       int OverlapLevel_in) :
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
  MpiRank_(0),
  NumMpiProcs_(1),
  NumMpiProcsPerSubdomain_(1),
  NumSubdomains_(1),
  SubdomainId_(0),
#endif
  IsInitialized_(false),
  IsComputed_(false),
  UseTranspose_(false),
  IsOverlapping_(false),
  OverlapLevel_(OverlapLevel_in),
  CombineMode_(Zero),
  Condest_(-1.0),
  ComputeCondest_(true),
  UseReordering_(false),
  ReorderingType_("none"),
  FilterSingletons_(false),
  NumInitialize_(0),
  NumCompute_(0),
  NumApplyInverse_(0),
  InitializeTime_(0.0),
  ComputeTime_(0.0),
  ApplyInverseTime_(0.0),
  InitializeFlops_(0.0),
  ComputeFlops_(0.0),
  ApplyInverseFlops_(0.0)
{
  // Construct a reference-counted pointer with the input matrix, don't manage the memory.
  Matrix_ = Teuchos::rcp( Matrix_in, false );

  if (Matrix_->Comm().NumProc() == 1)
    OverlapLevel_ = 0;

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

#ifdef IFPACK_NODE_AWARE_CODE
extern int ML_NODE_ID;
#endif

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

  Epetra_RowMatrix* MatrixPtr;

#ifndef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
# ifdef IFPACK_NODE_AWARE_CODE
/*
  sleep(3);
  if (Comm().MyPID() == 0) cout << "Printing out ovArowmap" << endl;
  Comm().Barrier();

  EpetraExt::BlockMapToMatrixMarketFile("ovArowmap",OverlappingMatrix_->RowMatrixRowMap());
  if (Comm().MyPID() == 0) cout << "Printing out ovAcolmap" << endl;
  Comm().Barrier();
  EpetraExt::BlockMapToMatrixMarketFile("ovAcolmap",OverlappingMatrix_->RowMatrixColMap());
  Comm().Barrier();
*/
/*
  EpetraExt::RowMatrixToMatlabFile("ovA",*OverlappingMatrix_);
  fprintf(stderr,"p %d n %d matrix file done\n",Comm().MyPID(),ML_NODE_ID);
  Comm().Barrier();
*/
  int nodeID;
  try{ nodeID = List_.get("ML node id",0);}
  catch(...){fprintf(stderr,"%s","Ifpack_AdditiveSchwarz<T>::Setup(): no parameter \"ML node id\"\n\n");
             cout << List_ << endl;}
# endif
#endif

  try{
  if (OverlappingMatrix_ != Teuchos::null)
  {
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
    if (NumMpiProcsPerSubdomain_ > 1) {
      LocalizedMatrix_ = Teuchos::rcp(new Ifpack_SubdomainFilter(OverlappingMatrix_, SubdomainId_));
    } else {
      LocalizedMatrix_ = Teuchos::rcp(new Ifpack_LocalFilter(OverlappingMatrix_));
    }
#else
#   ifdef IFPACK_NODE_AWARE_CODE
    Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(OverlappingMatrix_,nodeID); //FIXME
    LocalizedMatrix_ = Teuchos::rcp(tt);
    //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
#   else
    LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
#   endif
#endif
  }
  else
  {
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS

    if (NumMpiProcsPerSubdomain_ > 1) {
      LocalizedMatrix_ = Teuchos::rcp(new Ifpack_SubdomainFilter(Matrix_, SubdomainId_));
    } else {
      LocalizedMatrix_ = Teuchos::rcp(new Ifpack_LocalFilter(Matrix_));
    }
#else
#   ifdef IFPACK_NODE_AWARE_CODE
    Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(Matrix_,nodeID); //FIXME
    LocalizedMatrix_ = Teuchos::rcp(tt);
    //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
#   else
    LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
#   endif
#endif
  }
  }
  catch(...) {
     fprintf(stderr,"%s","AdditiveSchwarz Setup: problem creating local filter matrix.\n");
  }

  if (LocalizedMatrix_ == Teuchos::null)
    IFPACK_CHK_ERR(-5);

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

  if (UseReordering_) {

    // create reordering and compute it
    if (ReorderingType_ == "rcm")
      Reordering_ = Teuchos::rcp( new Ifpack_RCMReordering() );
    else if (ReorderingType_ == "metis")
      Reordering_ = Teuchos::rcp( new Ifpack_METISReordering() );
#ifdef HAVE_IFPACK_AMESOS   
    else if (ReorderingType_ == "amd" )
      Reordering_ = Teuchos::rcp( new Ifpack_AMDReordering() );
#endif
    else {
      cerr << "reordering type not correct (" << ReorderingType_ << ")" << endl;
      exit(EXIT_FAILURE);
    }
    if (Reordering_ == Teuchos::null) IFPACK_CHK_ERR(-5);

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

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

    if (ReorderedLocalizedMatrix_ == Teuchos::null) IFPACK_CHK_ERR(-5);

    MatrixPtr = &*ReorderedLocalizedMatrix_;
  }

#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
  // The subdomain matrix needs to be reindexed by Amesos so we need to make a CrsMatrix
  // and then reindex it with EpetraExt.
  // The reindexing is done here because this feature is only implemented in Amesos_Klu,
  // and it is desired to have reindexing with other direct solvers in the Amesos package

  SubdomainCrsMatrix_.reset(new Epetra_CrsMatrix(Copy, MatrixPtr->RowMatrixRowMap(), -1));
  Teuchos::RCP<Epetra_Import> tempImporter = Teuchos::rcp(new Epetra_Import(SubdomainCrsMatrix_->Map(), MatrixPtr->Map()));
  SubdomainCrsMatrix_->Import(*MatrixPtr, *tempImporter, Insert);
  SubdomainCrsMatrix_->FillComplete();

  if (NumMpiProcsPerSubdomain_ > 1) {
    tempMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->RowMap().NumGlobalElements(),
                                  SubdomainCrsMatrix_->RowMap().NumMyElements(),
                                  0, SubdomainCrsMatrix_->Comm()));
    tempRangeMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->OperatorRangeMap().NumGlobalElements(),
                                       SubdomainCrsMatrix_->OperatorRangeMap().NumMyElements(),
                                       0, SubdomainCrsMatrix_->Comm()));
    tempDomainMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->OperatorDomainMap().NumGlobalElements(),
                                        SubdomainCrsMatrix_->OperatorDomainMap().NumMyElements(),
                                        0, SubdomainCrsMatrix_->Comm()));

    SubdomainMatrixReindexer_.reset(new EpetraExt::CrsMatrix_Reindex(*tempMap_));
    DomainVectorReindexer_.reset(new EpetraExt::MultiVector_Reindex(*tempDomainMap_));
    RangeVectorReindexer_.reset(new EpetraExt::MultiVector_Reindex(*tempRangeMap_));

    ReindexedCrsMatrix_.reset(&((*SubdomainMatrixReindexer_)(*SubdomainCrsMatrix_)), false);

    MatrixPtr = &*ReindexedCrsMatrix_;

    Inverse_ = Teuchos::rcp( new T(&*ReindexedCrsMatrix_) );
  } else {
    Inverse_ = Teuchos::rcp( new T(&*SubdomainCrsMatrix_) );
  }
#else
  Inverse_ = Teuchos::rcp( new T(MatrixPtr) );
#endif

  if (Inverse_ == Teuchos::null)
    IFPACK_CHK_ERR(-5);

  return(0);
}

//==============================================================================
template<typename T>
int Ifpack_AdditiveSchwarz<T>::SetParameters(Teuchos::ParameterList& List_in)
{
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
  MpiRank_ = Matrix_->Comm().MyPID();
  NumMpiProcs_ = Matrix_->Comm().NumProc();
  NumMpiProcsPerSubdomain_ = List_in.get("subdomain: number-of-processors", 1);
  NumSubdomains_ = NumMpiProcs_ / NumMpiProcsPerSubdomain_;
  SubdomainId_ = MpiRank_ / NumMpiProcsPerSubdomain_;

  if (NumSubdomains_ == 1) {
      IsOverlapping_ = false;
  }

#endif
  // compute the condition number each time Compute() is invoked.
  ComputeCondest_ = List_in.get("schwarz: compute condest", ComputeCondest_);
  // combine mode
  if( Teuchos::ParameterEntry *combineModeEntry = List_in.getEntryPtr("schwarz: combine mode") )
  {
    if( typeid(std::string) == combineModeEntry->getAny().type() )
    {
      std::string mode = List_in.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
      {
        TEUCHOS_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_in,"schwarz: combine mode");
    }
  }
  else
  {
    // Make the default be a string to be consistent with the valid parameters!
    List_in.get("schwarz: combine mode","Zero");
  }
  // type of reordering
  ReorderingType_ = List_in.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_in.get("schwarz: filter singletons", FilterSingletons_);

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

  return(0);
}

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

  if (Time_ == Teuchos::null)
    Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );

  Time_->ResetStartTime();

  // compute the overlapping matrix if necessary
  if (IsOverlapping_) {
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
    if (NumMpiProcsPerSubdomain_ > 1) {
      OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, SubdomainId_) );
    } else {
      OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_));
    }
#else
#     ifdef IFPACK_NODE_AWARE_CODE
    int myNodeID;
    try{ myNodeID = List_.get("ML node id",-1);}
    catch(...){fprintf(stderr,"pid %d: no such entry (returned %d)\n",Comm().MyPID(),myNodeID);}
/*
      cout << "pid " << Comm().MyPID()
           << ": calling Ifpack_OverlappingRowMatrix with myNodeID = "
           << myNodeID << ", OverlapLevel_ = " << OverlapLevel_ << endl;
*/
    OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, myNodeID) );
#   else
    OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_) );
#   endif
#endif

    if (OverlappingMatrix_ == Teuchos::null) {
      IFPACK_CHK_ERR(-5);
    } 
  }

# ifdef IFPACK_NODE_AWARE_CODE
/*
  sleep(1);
  Comm().Barrier();
*/
# endif

  IFPACK_CHK_ERR(Setup());

# ifdef IFPACK_NODE_AWARE_CODE
/*
  sleep(1);
  Comm().Barrier();
*/
#endif

  if (Inverse_ == Teuchos::null)
    IFPACK_CHK_ERR(-5);

  if (LocalizedMatrix_ == Teuchos::null)
    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();

#ifdef IFPACK_FLOPCOUNTERS
  // 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;
#endif

  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();

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

  // 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_in)
{
  // store the flag -- it will be set in Initialize() if Inverse_ does not
  // exist.
  UseTranspose_ = UseTranspose_in;

  // If Inverse_ exists, pass it right now.
  if (Inverse_!=Teuchos::null)
    IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose_in));
  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();

  Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
  Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
  Teuchos::RefCountPtr<Epetra_MultiVector> Xtmp;

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

  // process overlap, may need to create vectors and import data
  if (IsOverlapping()) {
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
    if (OverlappingX == Teuchos::null) {
      OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(), X.NumVectors()) );
      if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
    } else assert(OverlappingX->NumVectors() == X.NumVectors());
    if (OverlappingY == Teuchos::null) {
      OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(), Y.NumVectors()) );
      if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
    } else assert(OverlappingY->NumVectors() == Y.NumVectors());
#else
#   ifdef IFPACK_NODE_AWARE_CODE
    if (OverlappingX == Teuchos::null) {
      OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
                                   X.NumVectors()) );
      if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
    } else assert(OverlappingX->NumVectors() == X.NumVectors());
    if (OverlappingY == Teuchos::null) {
      OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
                                     Y.NumVectors()) );
      if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
    } else assert(OverlappingY->NumVectors() == Y.NumVectors());
#else
    OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
                            X.NumVectors()) );
    OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
                            Y.NumVectors()) );
    if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
#   endif
#endif
    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 = Teuchos::rcp( new Epetra_MultiVector(X) );
    OverlappingX = Xtmp;
    OverlappingY = Teuchos::rcp( &Y, false );
  }

  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_) {
#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
        if (NumMpiProcsPerSubdomain_ > 1) {
          tempX_.reset(&((*RangeVectorReindexer_)(*OverlappingX)), false);
          tempY_.reset(&((*DomainVectorReindexer_)(*OverlappingY)), false);
          IFPACK_CHK_ERR(Inverse_->ApplyInverse(*tempX_,*tempY_));
        } else {
          IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX, *OverlappingY));
        }
#else
      IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX,*OverlappingY));
#endif
    }
    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_));
  }

#ifdef IFPACK_FLOPCOUNTERS
  // 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;
#endif

  // 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
{
#ifdef IFPACK_FLOPCOUNTERS
  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();
#endif

  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_->NumGlobalRows64() << endl;

#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
  os << endl;
  os << "================================================================================" << endl;
  os << "Subcommunicator stats" << endl;
  os << "Number of MPI processes in simulation: " << NumMpiProcs_ << endl;
  os << "Number of subdomains: " << NumSubdomains_ << endl;
  os << "Number of MPI processes per subdomain: " << NumMpiProcsPerSubdomain_ << endl;
#endif

  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() 
#ifdef IFPACK_FLOPCOUNTERS
     << "  " << std::setw(15) << 1.0e-6 * IF 
     << "  " << std::setw(15) << 1.0e-6 * IFT
#endif
     << endl;
  os << "Compute()       "   << std::setw(5) << NumCompute() 
     << "  " << std::setw(15) << ComputeTime()
#ifdef IFPACK_FLOPCOUNTERS
     << "  " << std::setw(15) << 1.0e-6 * CF
     << "  " << std::setw(15) << 1.0e-6 * CFT
#endif
     << endl;
  os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
     << "  " << std::setw(15) << ApplyInverseTime()
#ifdef IFPACK_FLOPCOUNTERS
     << "  " << std::setw(15) << 1.0e-6 * AF
     << "  " << std::setw(15) << 1.0e-6 * AFT
#endif
     << 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_in)
{
  if (!IsComputed()) // cannot compute right now
    return(-1.0);

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

  return(Condest_);
}

#endif // IFPACK_ADDITIVESCHWARZ_H