Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_IlukGraph.hpp
00001 /*@HEADER
00002  // ***********************************************************************
00003  // 
00004  //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005  //                 Copyright (2009) Sandia Corporation
00006  // 
00007  // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008  // license for use of this work by or on behalf of the U.S. Government.
00009  // 
00010  // This library is free software; you can redistribute it and/or modify
00011  // it under the terms of the GNU Lesser General Public License as
00012  // published by the Free Software Foundation; either version 2.1 of the
00013  // License, or (at your option) any later version.
00014  //  
00015  // This library is distributed in the hope that it will be useful, but
00016  // WITHOUT ANY WARRANTY; without even the implied warranty of
00017  // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018  // Lesser General Public License for more details.
00019  //  
00020  // You should have received a copy of the GNU Lesser General Public
00021  // License along with this library; if not, write to the Free Software
00022  // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023  // USA
00024  // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025  // 
00026  // ***********************************************************************
00027  //@HEADER
00028  */
00029 
00030 #ifndef IFPACK2_ILUK_GRAPH_HPP
00031 #define IFPACK2_ILUK_GRAPH_HPP
00032 
00033 #include <vector>
00034 #include <algorithm>
00035 
00036 #include <Ifpack2_ConfigDefs.hpp>
00037 #include <Teuchos_RCP.hpp>
00038 #include <Teuchos_ParameterList.hpp>
00039 #include <Teuchos_Describable.hpp>
00040 #include <Teuchos_CommHelpers.hpp>
00041 #include <Tpetra_CrsGraph.hpp>
00042 #include <Tpetra_Import.hpp>
00043 #include <Ifpack2_CreateOverlapGraph.hpp>
00044 #include <Ifpack2_Parameters.hpp>
00045 
00046 namespace Ifpack2 {
00047 
00048 
00050 
00069 template<class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
00070 class IlukGraph : public Teuchos::Describable
00071 {
00072 
00073 public:
00074   
00076 
00088   IlukGraph(const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > &Graph_in, int LevelFill_in, int LevelOverlap_in);
00089   
00091   IlukGraph(const IlukGraph<LocalOrdinal,GlobalOrdinal,Node> & Graph_in);
00092   
00094   virtual ~IlukGraph();
00095   
00097 
00102   void setParameters(const Teuchos::ParameterList& parameterlist);
00103   
00105 
00107   void constructFilledGraph();
00108   
00110 
00112   void constructOverlapGraph();
00113   
00115   int getLevelFill() const {return(LevelFill_);}
00116   
00118   int getLevelOverlap() const {return(LevelOverlap_);}
00119   
00121   const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& getL_Graph() const {return(L_Graph_);}
00122   
00124   const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& getU_Graph() const {return(U_Graph_);}
00125   
00127   const Teuchos::RCP<Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& getOverlapGraph() const  {return(OverlapGraph_);}
00128 
00130   size_t getNumGlobalDiagonals() const { return NumGlobalDiagonals_; }
00131 
00132 private:
00133   typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> TpetraMapType;
00134   typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> TpetraCrsGraphType;
00135   
00136   Teuchos::RCP<const TpetraCrsGraphType > Graph_;
00137   Teuchos::RCP<const TpetraCrsGraphType > OverlapGraph_;
00138   int LevelFill_;
00139   int LevelOverlap_;
00140   Teuchos::RCP<TpetraCrsGraphType > L_Graph_;
00141   Teuchos::RCP<TpetraCrsGraphType > U_Graph_;
00142   size_t NumMyDiagonals_;
00143   size_t NumGlobalDiagonals_;
00144 };
00145 
00146 //==============================================================================
00147 template<class LocalOrdinal, class GlobalOrdinal, class Node>
00148 IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::IlukGraph(const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& Graph_in, int LevelFill_in, int LevelOverlap_in)
00149 : Graph_(Graph_in),
00150   OverlapGraph_(),
00151   LevelFill_(LevelFill_in),
00152   LevelOverlap_(LevelOverlap_in),
00153   NumMyDiagonals_(0),
00154   NumGlobalDiagonals_(0)
00155 {
00156 }
00157 
00158 //==============================================================================
00159 template<class LocalOrdinal, class GlobalOrdinal, class Node>
00160 IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::IlukGraph(const IlukGraph<LocalOrdinal,GlobalOrdinal,Node> & Graph_in) 
00161 : Graph_(Graph_in.Graph_),
00162   OverlapGraph_(Graph_in.OverlapGraph_),
00163   LevelFill_(Graph_in.LevelFill_),
00164   LevelOverlap_(Graph_in.LevelOverlap_),
00165   L_Graph_(),
00166   U_Graph_(),
00167   NumMyDiagonals_(0),
00168   NumGlobalDiagonals_(0)
00169 {
00170   TpetraCrsGraphType & L_Graph_In = Graph_in.L_Graph();
00171   TpetraCrsGraphType & U_Graph_In = Graph_in.U_Graph();
00172   L_Graph_ = Teuchos::rcp( new TpetraCrsGraphType(L_Graph_In) );
00173   U_Graph_ = Teuchos::rcp( new TpetraCrsGraphType(U_Graph_In) );
00174 }
00175 
00176 //==============================================================================
00177 template<class LocalOrdinal, class GlobalOrdinal, class Node>
00178 IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::~IlukGraph()
00179 {
00180 }
00181 
00182 //==============================================================================
00183 template<class LocalOrdinal, class GlobalOrdinal, class Node>
00184 void IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::setParameters(const Teuchos::ParameterList& parameterlist)
00185 {
00186   getParameter(parameterlist, "iluk level-of-fill", LevelFill_);
00187   getParameter(parameterlist, "iluk level-of-overlap", LevelOverlap_);
00188 }
00189 
00190 //==============================================================================
00191 template<class LocalOrdinal, class GlobalOrdinal, class Node>
00192 void IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::constructOverlapGraph() {
00193   
00194   if (OverlapGraph_ == Teuchos::null) {
00195     OverlapGraph_ = CreateOverlapGraph(Graph_, LevelOverlap_);
00196   }
00197 }
00198 
00199 //==============================================================================
00200 template<class LocalOrdinal, class GlobalOrdinal, class Node>
00201 void IlukGraph<LocalOrdinal,GlobalOrdinal,Node>::constructFilledGraph() {
00202   size_t NumIn, NumL, NumU;
00203   bool DiagFound;
00204  
00205   constructOverlapGraph();
00206  
00207   L_Graph_ = Teuchos::rcp( new TpetraCrsGraphType(OverlapGraph_->getRowMap(), OverlapGraph_->getRowMap(),  0));
00208   U_Graph_ = Teuchos::rcp( new TpetraCrsGraphType(OverlapGraph_->getRowMap(), OverlapGraph_->getRowMap(),  0));
00209  
00210  
00211   // Get Maximum Row length
00212   int MaxNumIndices = OverlapGraph_->getNodeMaxNumRowEntries();
00213  
00214   Teuchos::Array<LocalOrdinal> L(MaxNumIndices);
00215   Teuchos::Array<LocalOrdinal> U(MaxNumIndices);
00216  
00217   // First we copy the user's graph into L and U, regardless of fill level
00218  
00219   int NumMyRows = OverlapGraph_->getRowMap()->getNodeNumElements();
00220   NumMyDiagonals_ = 0;
00221 
00222   for (int i=0; i< NumMyRows; i++) {
00223  
00224     Teuchos::ArrayView<const LocalOrdinal> my_indices;
00225     OverlapGraph_->getLocalRowView(i,my_indices);
00226  
00227     // Split into L and U (we don't assume that indices are ordered).
00228     
00229     NumL = 0; 
00230     NumU = 0; 
00231     DiagFound = false;
00232     NumIn = my_indices.size();
00233 
00234     for (size_t j=0; j< NumIn; j++) {
00235       LocalOrdinal k = my_indices[j];
00236       
00237       if (k<NumMyRows) { // Ignore column elements that are not in the square matrix
00238         
00239         if (k==i) DiagFound = true;
00240         
00241         else if (k < i) {
00242           L[NumL] = k;
00243           NumL++;
00244         }
00245         else {
00246           U[NumU] = k;
00247           NumU++;
00248         }
00249       }
00250     }
00251 
00252     // Check in things for this row of L and U
00253  
00254     if (DiagFound) ++NumMyDiagonals_;
00255     if (NumL) {
00256       Teuchos::ArrayView<LocalOrdinal> Lview(&L[0], NumL);
00257       L_Graph_->insertLocalIndices(i, Lview );
00258     }
00259     if (NumU) {
00260       Teuchos::ArrayView<LocalOrdinal> Uview(&U[0], NumU);
00261       U_Graph_->insertLocalIndices(i, Uview );
00262     }
00263   }
00264  
00265   if (LevelFill_ > 0) {
00266     
00267     // Complete Fill steps
00268     Teuchos::RCP<const TpetraMapType> L_DomainMap = OverlapGraph_->getRowMap();
00269     Teuchos::RCP<const TpetraMapType> L_RangeMap = Graph_->getRangeMap();
00270     Teuchos::RCP<const TpetraMapType> U_DomainMap = Graph_->getDomainMap();
00271     Teuchos::RCP<const TpetraMapType> U_RangeMap = OverlapGraph_->getRowMap();
00272     Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(new Teuchos::ParameterList());
00273     params->set("Optimize Storage",false);
00274     L_Graph_->fillComplete(L_DomainMap, L_RangeMap, params);
00275     U_Graph_->fillComplete(U_DomainMap, U_RangeMap, params);
00276     L_Graph_->resumeFill();
00277     U_Graph_->resumeFill();
00278     
00279     // At this point L_Graph and U_Graph are filled with the pattern of input graph, 
00280     // sorted and have redundant indices (if any) removed.  Indices are zero based.
00281     // LevelFill is greater than zero, so continue...
00282     
00283     int MaxRC = NumMyRows;
00284     std::vector<std::vector<int> > Levels(MaxRC);
00285     std::vector<int> LinkList(MaxRC);
00286     std::vector<int> CurrentLevel(MaxRC);
00287     Teuchos::Array<LocalOrdinal> CurrentRow(MaxRC+1);
00288     std::vector<int> LevelsRowU(MaxRC);
00289  
00290     for (int i=0; i<NumMyRows; i++)
00291     {
00292       int First, Next;
00293 
00294       // copy column indices of row into workspace and sort them
00295 
00296       size_t LenL = L_Graph_->getNumEntriesInLocalRow(i);
00297       size_t LenU = U_Graph_->getNumEntriesInLocalRow(i);
00298       size_t Len = LenL + LenU + 1;
00299 
00300       CurrentRow.resize(Len);
00301 
00302       L_Graph_->getLocalRowCopy(i, CurrentRow(), LenL);      // Get L Indices
00303       CurrentRow[LenL] = i;                                     // Put in Diagonal
00304       if (LenU > 0) {
00305         Teuchos::ArrayView<LocalOrdinal> URowView(&CurrentRow[LenL+1], LenU);
00306         // Get U Indices
00307         U_Graph_->getLocalRowCopy(i, URowView, LenU);
00308       }
00309 
00310       // Construct linked list for current row
00311       
00312       for (size_t j=0; j<Len-1; j++) {
00313         LinkList[CurrentRow[j]] = CurrentRow[j+1];
00314         CurrentLevel[CurrentRow[j]] = 0;
00315       }
00316       
00317       LinkList[CurrentRow[Len-1]] = NumMyRows;
00318       CurrentLevel[CurrentRow[Len-1]] = 0;
00319       
00320       // Merge List with rows in U
00321       
00322       First = CurrentRow[0];
00323       Next = First;
00324       while (Next < i)
00325       {
00326         int PrevInList = Next;
00327         int NextInList = LinkList[Next];
00328         int RowU = Next;
00329         // Get Indices for this row of U
00330         Teuchos::ArrayView<const LocalOrdinal> IndicesU;
00331         U_Graph_->getLocalRowView(RowU,IndicesU);
00332         int LengthRowU = IndicesU.size();
00333         
00334         int ii;
00335         
00336         // Scan RowU
00337         
00338         for (ii=0; ii<LengthRowU; /*nop*/)
00339         {
00340           int CurInList = IndicesU[ii];
00341           if (CurInList < NextInList)
00342           {
00343             // new fill-in
00344             int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00345             if (NewLevel <= LevelFill_)
00346             {
00347               LinkList[PrevInList]  = CurInList;
00348               LinkList[CurInList] = NextInList;
00349               PrevInList = CurInList;
00350               CurrentLevel[CurInList] = NewLevel;
00351             }
00352             ii++;
00353           }
00354           else if (CurInList == NextInList)
00355           {
00356             PrevInList = NextInList;
00357             NextInList = LinkList[PrevInList];
00358             int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00359             CurrentLevel[CurInList] = std::min(CurrentLevel[CurInList], NewLevel);
00360             ii++;
00361           }
00362           else // (CurInList > NextInList)
00363           {
00364             PrevInList = NextInList;
00365             NextInList = LinkList[PrevInList];
00366           }
00367         }
00368         Next = LinkList[Next];
00369       }
00370       
00371       // Put pattern into L and U
00372       
00373       CurrentRow.resize(0);
00374       
00375       Next = First;
00376       
00377       // Lower
00378       
00379       while (Next < i) {    
00380         CurrentRow.push_back(Next);
00381         Next = LinkList[Next];
00382       }
00383 
00384       L_Graph_->removeLocalIndices(i); // Delete current set of Indices
00385       if (CurrentRow.size() > 0) {
00386         L_Graph_->insertLocalIndices(i, CurrentRow());
00387       }
00388  
00389       // Diagonal
00390       
00391       TEUCHOS_TEST_FOR_EXCEPTION(Next != i, std::runtime_error,
00392                          "Ifpack2::IlukGraph::constructFilledGraph: FATAL: U has zero diagonal")
00393 
00394       LevelsRowU[0] = CurrentLevel[Next];
00395       Next = LinkList[Next];
00396       
00397       // Upper
00398       
00399       CurrentRow.resize(0);
00400       LenU = 0;
00401       
00402       while (Next < NumMyRows) {
00403         LevelsRowU[LenU+1] = CurrentLevel[Next];
00404         CurrentRow.push_back(Next);
00405         ++LenU;
00406         Next = LinkList[Next];
00407       }
00408       
00409       U_Graph_->removeLocalIndices(i); // Delete current set of Indices
00410       if (LenU > 0) {
00411         U_Graph_->insertLocalIndices(i, CurrentRow());
00412       }
00413  
00414       // Allocate and fill Level info for this row
00415       Levels[i] = std::vector<int>(LenU+1);
00416       for (size_t jj=0; jj<LenU+1; jj++) {
00417         Levels[i][jj] = LevelsRowU[jj];
00418       }
00419     }
00420   }    
00421   
00422   // Complete Fill steps
00423   Teuchos::RCP<const TpetraMapType> L_DomainMap = OverlapGraph_->getRowMap();
00424   Teuchos::RCP<const TpetraMapType> L_RangeMap  = Graph_->getRangeMap();
00425   Teuchos::RCP<const TpetraMapType> U_DomainMap = Graph_->getDomainMap();
00426   Teuchos::RCP<const TpetraMapType> U_RangeMap  = OverlapGraph_->getRowMap();
00427   L_Graph_->fillComplete(L_DomainMap, L_RangeMap);//DoOptimizeStorage is default here...
00428   U_Graph_->fillComplete(U_DomainMap, U_RangeMap);//DoOptimizeStorage is default here...
00429 
00430   Teuchos::reduceAll<int,size_t>(*(L_DomainMap->getComm()), Teuchos::REDUCE_SUM, 1, &NumMyDiagonals_, &NumGlobalDiagonals_);
00431 }
00432 
00433 } // namespace Ifpack2
00434 
00435 #endif /* IFPACK2_ILUK_GRAPH_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends