Ifpack2_IlukGraph.hpp

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

#include <vector>
#include <algorithm>

#include <Ifpack2_ConfigDefs.hpp>
#include <Teuchos_RCP.hpp>
#include <Teuchos_ParameterList.hpp>
#include <Teuchos_Describable.hpp>
#include <Teuchos_CommHelpers.hpp>
#include <Tpetra_CrsGraph.hpp>
#include <Tpetra_Import.hpp>
#include <Ifpack2_CreateOverlapGraph.hpp>
#include <Ifpack2_Parameters.hpp>

namespace Ifpack2 {



template<class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
class IlukGraph : public Teuchos::Describable
{

public:
  

  IlukGraph(const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > &Graph_in, int LevelFill_in, int LevelOverlap_in);
  
  IlukGraph(const IlukGraph<LocalOrdinal,GlobalOrdinal,Node> & Graph_in);
  
  virtual ~IlukGraph();
  

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

  void constructFilledGraph();
  

  void constructOverlapGraph();
  
  int getLevelFill() const {return(LevelFill_);}
  
  int getLevelOverlap() const {return(LevelOverlap_);}
  
  const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& getL_Graph() const {return(L_Graph_);}
  
  const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& getU_Graph() const {return(U_Graph_);}
  
  const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& getOverlapGraph() const  {return(OverlapGraph_);}

  size_t getNumGlobalDiagonals() const { return NumGlobalDiagonals_; }

private:
  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpetraMapType;
  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> TpetraCrsGraphType;
  
  Teuchos::RCP<const TpetraCrsGraphType > Graph_;
  Teuchos::RCP<const TpetraCrsGraphType > OverlapGraph_;
  int LevelFill_;
  int LevelOverlap_;
  Teuchos::RCP<TpetraCrsGraphType > L_Graph_;
  Teuchos::RCP<TpetraCrsGraphType > U_Graph_;
  size_t NumMyDiagonals_;
  size_t NumGlobalDiagonals_;
};

//==============================================================================
template<class LocalOrdinal, class GlobalOrdinal, class Node>
IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::IlukGraph(const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& Graph_in, int LevelFill_in, int LevelOverlap_in)
: Graph_(Graph_in),
  OverlapGraph_(),
  LevelFill_(LevelFill_in),
  LevelOverlap_(LevelOverlap_in),
  NumMyDiagonals_(0),
  NumGlobalDiagonals_(0)
{
}

//==============================================================================
template<class LocalOrdinal, class GlobalOrdinal, class Node>
IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::IlukGraph(const IlukGraph<LocalOrdinal,GlobalOrdinal,Node> & Graph_in) 
: Graph_(Graph_in.Graph_),
  OverlapGraph_(Graph_in.OverlapGraph_),
  LevelFill_(Graph_in.LevelFill_),
  LevelOverlap_(Graph_in.LevelOverlap_),
  L_Graph_(),
  U_Graph_(),
  NumMyDiagonals_(0),
  NumGlobalDiagonals_(0)
{
  TpetraCrsGraphType & L_Graph_In = Graph_in.L_Graph();
  TpetraCrsGraphType & U_Graph_In = Graph_in.U_Graph();
  L_Graph_ = Teuchos::rcp( new TpetraCrsGraphType(L_Graph_In) );
  U_Graph_ = Teuchos::rcp( new TpetraCrsGraphType(U_Graph_In) );
}

//==============================================================================
template<class LocalOrdinal, class GlobalOrdinal, class Node>
IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::~IlukGraph()
{
}

//==============================================================================
template<class LocalOrdinal, class GlobalOrdinal, class Node>
void IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::setParameters(const Teuchos::ParameterList& parameterlist)
{
  getParameter(parameterlist, "iluk level-of-fill", LevelFill_);
  getParameter(parameterlist, "iluk level-of-overlap", LevelOverlap_);
}

//==============================================================================
template<class LocalOrdinal, class GlobalOrdinal, class Node>
void IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::constructOverlapGraph() {
  
  if (OverlapGraph_ == Teuchos::null) {
    OverlapGraph_ = CreateOverlapGraph(Graph_, LevelOverlap_);
  }
}

//==============================================================================
template<class LocalOrdinal, class GlobalOrdinal, class Node>
void IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::constructFilledGraph() {
  size_t NumIn, NumL, NumU;
  bool DiagFound;
 
  constructOverlapGraph();
 
  L_Graph_ = Teuchos::rcp( new TpetraCrsGraphType(OverlapGraph_->getRowMap(), OverlapGraph_->getRowMap(),  0));
  U_Graph_ = Teuchos::rcp( new TpetraCrsGraphType(OverlapGraph_->getRowMap(), OverlapGraph_->getRowMap(),  0));
 
 
  // Get Maximum Row length
  int MaxNumIndices = OverlapGraph_->getNodeMaxNumRowEntries();
 
  Teuchos::Array<LocalOrdinal> L(MaxNumIndices);
  Teuchos::Array<LocalOrdinal> U(MaxNumIndices);
 
  // First we copy the user's graph into L and U, regardless of fill level
 
  int NumMyRows = OverlapGraph_->getRowMap()->getNodeNumElements();
  NumMyDiagonals_ = 0;

  for (int i=0; i< NumMyRows; i++) {
 
    Teuchos::ArrayView<const LocalOrdinal> my_indices;
    OverlapGraph_->getLocalRowView(i,my_indices);
 
    // Split into L and U (we don't assume that indices are ordered).
    
    NumL = 0; 
    NumU = 0; 
    DiagFound = false;
    NumIn = my_indices.size();

    for (size_t j=0; j< NumIn; j++) {
      LocalOrdinal k = my_indices[j];
      
      if (k<NumMyRows) { // Ignore column elements that are not in the square matrix
        
        if (k==i) DiagFound = true;
        
        else if (k < i) {
          L[NumL] = k;
          NumL++;
        }
        else {
          U[NumU] = k;
          NumU++;
        }
      }
    }

    // Check in things for this row of L and U
 
    if (DiagFound) ++NumMyDiagonals_;
    if (NumL) {
      Teuchos::ArrayView<LocalOrdinal> Lview(&L[0], NumL);
      L_Graph_->insertLocalIndices(i, Lview );
    }
    if (NumU) {
      Teuchos::ArrayView<LocalOrdinal> Uview(&U[0], NumU);
      U_Graph_->insertLocalIndices(i, Uview );
    }
  }
 
  if (LevelFill_ > 0) {
    
    // Complete Fill steps
    Teuchos::RCP<const TpetraMapType> L_DomainMap = OverlapGraph_->getRowMap();
    Teuchos::RCP<const TpetraMapType> L_RangeMap = Graph_->getRangeMap();
    Teuchos::RCP<const TpetraMapType> U_DomainMap = Graph_->getDomainMap();
    Teuchos::RCP<const TpetraMapType> U_RangeMap = OverlapGraph_->getRowMap();
    Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(new Teuchos::ParameterList());
    params->set("Optimize Storage",false);
    L_Graph_->fillComplete(L_DomainMap, L_RangeMap, params);
    U_Graph_->fillComplete(U_DomainMap, U_RangeMap, params);
    L_Graph_->resumeFill();
    U_Graph_->resumeFill();
    
    // At this point L_Graph and U_Graph are filled with the pattern of input graph, 
    // sorted and have redundant indices (if any) removed.  Indices are zero based.
    // LevelFill is greater than zero, so continue...
    
    int MaxRC = NumMyRows;
    std::vector<std::vector<int> > Levels(MaxRC);
    std::vector<int> LinkList(MaxRC);
    std::vector<int> CurrentLevel(MaxRC);
    Teuchos::Array<LocalOrdinal> CurrentRow(MaxRC+1);
    std::vector<int> LevelsRowU(MaxRC);
 
    for (int i=0; i<NumMyRows; i++)
    {
      int First, Next;

      // copy column indices of row into workspace and sort them

      size_t LenL = L_Graph_->getNumEntriesInLocalRow(i);
      size_t LenU = U_Graph_->getNumEntriesInLocalRow(i);
      size_t Len = LenL + LenU + 1;

      CurrentRow.resize(Len);

      L_Graph_->getLocalRowCopy(i, CurrentRow(), LenL);      // Get L Indices
      CurrentRow[LenL] = i;                                     // Put in Diagonal
      if (LenU > 0) {
        Teuchos::ArrayView<LocalOrdinal> URowView(&CurrentRow[LenL+1], LenU);
        // Get U Indices
        U_Graph_->getLocalRowCopy(i, URowView, LenU);
      }

      // Construct linked list for current row
      
      for (size_t j=0; j<Len-1; j++) {
        LinkList[CurrentRow[j]] = CurrentRow[j+1];
        CurrentLevel[CurrentRow[j]] = 0;
      }
      
      LinkList[CurrentRow[Len-1]] = NumMyRows;
      CurrentLevel[CurrentRow[Len-1]] = 0;
      
      // Merge List with rows in U
      
      First = CurrentRow[0];
      Next = First;
      while (Next < i)
      {
        int PrevInList = Next;
        int NextInList = LinkList[Next];
        int RowU = Next;
        // Get Indices for this row of U
        Teuchos::ArrayView<const LocalOrdinal> IndicesU;
        U_Graph_->getLocalRowView(RowU,IndicesU);
        int LengthRowU = IndicesU.size();
        
        int ii;
        
        // Scan RowU
        
        for (ii=0; ii<LengthRowU; /*nop*/)
        {
          int CurInList = IndicesU[ii];
          if (CurInList < NextInList)
          {
            // new fill-in
            int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
            if (NewLevel <= LevelFill_)
            {
              LinkList[PrevInList]  = CurInList;
              LinkList[CurInList] = NextInList;
              PrevInList = CurInList;
              CurrentLevel[CurInList] = NewLevel;
            }
            ii++;
          }
          else if (CurInList == NextInList)
          {
            PrevInList = NextInList;
            NextInList = LinkList[PrevInList];
            int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
            CurrentLevel[CurInList] = std::min(CurrentLevel[CurInList], NewLevel);
            ii++;
          }
          else // (CurInList > NextInList)
          {
            PrevInList = NextInList;
            NextInList = LinkList[PrevInList];
          }
        }
        Next = LinkList[Next];
      }
      
      // Put pattern into L and U
      
      CurrentRow.resize(0);
      
      Next = First;
      
      // Lower
      
      while (Next < i) {    
        CurrentRow.push_back(Next);
        Next = LinkList[Next];
      }

      L_Graph_->removeLocalIndices(i); // Delete current set of Indices
      if (CurrentRow.size() > 0) {
        L_Graph_->insertLocalIndices(i, CurrentRow());
      }
 
      // Diagonal
      
      TEUCHOS_TEST_FOR_EXCEPTION(Next != i, std::runtime_error,
                         "Ifpack2::IlukGraph::constructFilledGraph: FATAL: U has zero diagonal")

      LevelsRowU[0] = CurrentLevel[Next];
      Next = LinkList[Next];
      
      // Upper
      
      CurrentRow.resize(0);
      LenU = 0;
      
      while (Next < NumMyRows) {
        LevelsRowU[LenU+1] = CurrentLevel[Next];
        CurrentRow.push_back(Next);
        ++LenU;
        Next = LinkList[Next];
      }
      
      U_Graph_->removeLocalIndices(i); // Delete current set of Indices
      if (LenU > 0) {
        U_Graph_->insertLocalIndices(i, CurrentRow());
      }
 
      // Allocate and fill Level info for this row
      Levels[i] = std::vector<int>(LenU+1);
      for (size_t jj=0; jj<LenU+1; jj++) {
        Levels[i][jj] = LevelsRowU[jj];
      }
    }
  }    
  
  // Complete Fill steps
  Teuchos::RCP<const TpetraMapType> L_DomainMap = OverlapGraph_->getRowMap();
  Teuchos::RCP<const TpetraMapType> L_RangeMap  = Graph_->getRangeMap();
  Teuchos::RCP<const TpetraMapType> U_DomainMap = Graph_->getDomainMap();
  Teuchos::RCP<const TpetraMapType> U_RangeMap  = OverlapGraph_->getRowMap();
  L_Graph_->fillComplete(L_DomainMap, L_RangeMap);//DoOptimizeStorage is default here...
  U_Graph_->fillComplete(U_DomainMap, U_RangeMap);//DoOptimizeStorage is default here...

  Teuchos::reduceAll<int,size_t>(*(L_DomainMap->getComm()), Teuchos::REDUCE_SUM, 1, &NumMyDiagonals_, &NumGlobalDiagonals_);
}

} // namespace Ifpack2

#endif /* IFPACK2_ILUK_GRAPH_HPP */