Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_IlukGraph.hpp
Go to the documentation of this file.
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041  */
00042 
00048 
00049 #ifndef IFPACK2_ILUK_GRAPH_HPP
00050 #define IFPACK2_ILUK_GRAPH_HPP
00051 
00052 #include <algorithm>
00053 #include <vector>
00054 
00055 #include <Ifpack2_ConfigDefs.hpp>
00056 #include <Teuchos_ParameterList.hpp>
00057 #include <Teuchos_CommHelpers.hpp>
00058 #include <Tpetra_CrsGraph.hpp>
00059 #include <Tpetra_Import.hpp>
00060 #include <Ifpack2_CreateOverlapGraph.hpp>
00061 #include <Ifpack2_Parameters.hpp>
00062 
00063 namespace Ifpack2 {
00064 
00096 template<class GraphType>
00097 class IlukGraph : public Teuchos::Describable {
00098 public:
00099   typedef typename GraphType::local_ordinal_type local_ordinal_type;
00100   typedef typename GraphType::global_ordinal_type global_ordinal_type;
00101   typedef typename GraphType::node_type node_type;
00102 
00104   typedef Tpetra::RowGraph<local_ordinal_type,
00105                            global_ordinal_type,
00106                            node_type> row_graph_type;
00108   typedef Tpetra::CrsGraph<local_ordinal_type,
00109                            global_ordinal_type,
00110                            node_type> crs_graph_type;
00111 
00122   IlukGraph (const Teuchos::RCP<const GraphType>& G,
00123              const int levelFill,
00124              const int levelOverlap);
00125 
00127   virtual ~IlukGraph ();
00128 
00134   void setParameters (const Teuchos::ParameterList& parameterlist);
00135 
00147   void initialize ();
00148 
00150   int getLevelFill () const { return LevelFill_; }
00151 
00153   int getLevelOverlap () const { return LevelOverlap_; }
00154 
00156   Teuchos::RCP<crs_graph_type> getL_Graph () const {
00157     return L_Graph_;
00158   }
00159 
00161   Teuchos::RCP<crs_graph_type> getU_Graph () const {
00162     return U_Graph_;
00163   }
00164 
00166   Teuchos::RCP<crs_graph_type> getOverlapGraph () const {
00167     return OverlapGraph_;
00168   }
00169 
00171   size_t getNumGlobalDiagonals() const { return NumGlobalDiagonals_; }
00172 
00173 private:
00174   typedef typename GraphType::map_type map_type;
00175 
00184   IlukGraph (const IlukGraph<GraphType>&);
00185 
00194   IlukGraph& operator= (const IlukGraph<GraphType>&);
00195 
00197   void constructOverlapGraph();
00198 
00199   Teuchos::RCP<const GraphType> Graph_;
00200   Teuchos::RCP<const crs_graph_type> OverlapGraph_;
00201   int LevelFill_;
00202   int LevelOverlap_;
00203   Teuchos::RCP<crs_graph_type> L_Graph_;
00204   Teuchos::RCP<crs_graph_type> U_Graph_;
00205   size_t NumMyDiagonals_;
00206   size_t NumGlobalDiagonals_;
00207 };
00208 
00209 
00210 template<class GraphType>
00211 IlukGraph<GraphType>::
00212 IlukGraph (const Teuchos::RCP<const GraphType>& G,
00213            const int levelFill,
00214            const int levelOverlap)
00215   : Graph_ (G),
00216     LevelFill_ (levelFill),
00217     LevelOverlap_ (levelOverlap),
00218     NumMyDiagonals_ (0),
00219     NumGlobalDiagonals_ (0)
00220 {}
00221 
00222 
00223 template<class GraphType>
00224 IlukGraph<GraphType>::~IlukGraph()
00225 {}
00226 
00227 
00228 template<class GraphType>
00229 void IlukGraph<GraphType>::
00230 setParameters (const Teuchos::ParameterList& parameterlist)
00231 {
00232   getParameter (parameterlist, "iluk level-of-fill", LevelFill_);
00233   getParameter (parameterlist, "iluk level-of-overlap", LevelOverlap_);
00234 }
00235 
00236 
00237 template<class GraphType>
00238 void IlukGraph<GraphType>::constructOverlapGraph () {
00239   // FIXME (mfh 22 Dec 2013) This won't do if we want
00240   // RILUK::initialize() to do the right thing (that is,
00241   // unconditionally recompute the "symbolic factorization").
00242   if (OverlapGraph_ == Teuchos::null) {
00243     OverlapGraph_ = createOverlapGraph<GraphType> (Graph_, LevelOverlap_);
00244   }
00245 }
00246 
00247 
00248 template<class GraphType>
00249 void IlukGraph<GraphType>::initialize()
00250 {
00251   using Teuchos::Array;
00252   using Teuchos::ArrayView;
00253   using Teuchos::RCP;
00254   using Teuchos::rcp;
00255   using Teuchos::REDUCE_SUM;
00256   using Teuchos::reduceAll;
00257 
00258   size_t NumIn, NumL, NumU;
00259   bool DiagFound;
00260 
00261   constructOverlapGraph();
00262 
00263   L_Graph_ = rcp (new crs_graph_type (OverlapGraph_->getRowMap (),
00264                                       OverlapGraph_->getRowMap (), 0));
00265   U_Graph_ = rcp (new crs_graph_type (OverlapGraph_->getRowMap (),
00266                                       OverlapGraph_->getRowMap (), 0));
00267 
00268   // Get Maximum Row length
00269   const int MaxNumIndices = OverlapGraph_->getNodeMaxNumRowEntries ();
00270 
00271   Array<local_ordinal_type> L (MaxNumIndices);
00272   Array<local_ordinal_type> U (MaxNumIndices);
00273 
00274   // First we copy the user's graph into L and U, regardless of fill level
00275 
00276   // FIXME (mfh 23 Dec 2013) Use size_t or whatever
00277   // getNodeNumElements() returns, instead of ptrdiff_t.
00278   const int NumMyRows = OverlapGraph_->getRowMap ()->getNodeNumElements ();
00279   NumMyDiagonals_ = 0;
00280 
00281   for (int i = 0; i< NumMyRows; ++i) {
00282     ArrayView<const local_ordinal_type> my_indices;
00283     OverlapGraph_->getLocalRowView (i, my_indices);
00284 
00285     // Split into L and U (we don't assume that indices are ordered).
00286 
00287     NumL = 0;
00288     NumU = 0;
00289     DiagFound = false;
00290     NumIn = my_indices.size();
00291 
00292     for (size_t j = 0; j < NumIn; ++j) {
00293       const local_ordinal_type k = my_indices[j];
00294 
00295       if (k<NumMyRows) { // Ignore column elements that are not in the square matrix
00296 
00297         if (k==i) {
00298           DiagFound = true;
00299         }
00300         else if (k < i) {
00301           L[NumL] = k;
00302           NumL++;
00303         }
00304         else {
00305           U[NumU] = k;
00306           NumU++;
00307         }
00308       }
00309     }
00310 
00311     // Check in things for this row of L and U
00312 
00313     if (DiagFound) {
00314       ++NumMyDiagonals_;
00315     }
00316     if (NumL) {
00317       ArrayView<const local_ordinal_type> L_view = L.view (0, NumL);
00318       L_Graph_->insertLocalIndices (i, L_view);
00319     }
00320     if (NumU) {
00321       ArrayView<const local_ordinal_type> U_view = U.view (0, NumU);
00322       U_Graph_->insertLocalIndices (i, U_view);
00323     }
00324   }
00325 
00326   if (LevelFill_ > 0) {
00327     // Complete Fill steps
00328     RCP<const map_type> L_DomainMap = OverlapGraph_->getRowMap ();
00329     RCP<const map_type> L_RangeMap = Graph_->getRangeMap ();
00330     RCP<const map_type> U_DomainMap = Graph_->getDomainMap ();
00331     RCP<const map_type> U_RangeMap = OverlapGraph_->getRowMap ();
00332     RCP<Teuchos::ParameterList> params = Teuchos::parameterList ();
00333     params->set ("Optimize Storage",false);
00334     L_Graph_->fillComplete (L_DomainMap, L_RangeMap, params);
00335     U_Graph_->fillComplete (U_DomainMap, U_RangeMap, params);
00336     L_Graph_->resumeFill ();
00337     U_Graph_->resumeFill ();
00338 
00339     // At this point L_Graph and U_Graph are filled with the pattern of input graph,
00340     // sorted and have redundant indices (if any) removed.  Indices are zero based.
00341     // LevelFill is greater than zero, so continue...
00342 
00343     int MaxRC = NumMyRows;
00344     std::vector<std::vector<int> > Levels(MaxRC);
00345     std::vector<int> LinkList(MaxRC);
00346     std::vector<int> CurrentLevel(MaxRC);
00347     Array<local_ordinal_type> CurrentRow (MaxRC + 1);
00348     std::vector<int> LevelsRowU(MaxRC);
00349 
00350     for (int i = 0; i < NumMyRows; ++i) {
00351       int First, Next;
00352 
00353       // copy column indices of row into workspace and sort them
00354 
00355       size_t LenL = L_Graph_->getNumEntriesInLocalRow(i);
00356       size_t LenU = U_Graph_->getNumEntriesInLocalRow(i);
00357       size_t Len = LenL + LenU + 1;
00358 
00359       CurrentRow.resize(Len);
00360 
00361       L_Graph_->getLocalRowCopy(i, CurrentRow(), LenL);      // Get L Indices
00362       CurrentRow[LenL] = i;                                     // Put in Diagonal
00363       if (LenU > 0) {
00364         ArrayView<local_ordinal_type> URowView = CurrentRow.view (LenL+1, LenU);
00365         // Get U Indices
00366         U_Graph_->getLocalRowCopy (i, URowView, LenU);
00367       }
00368 
00369       // Construct linked list for current row
00370 
00371       for (size_t j=0; j<Len-1; j++) {
00372         LinkList[CurrentRow[j]] = CurrentRow[j+1];
00373         CurrentLevel[CurrentRow[j]] = 0;
00374       }
00375 
00376       LinkList[CurrentRow[Len-1]] = NumMyRows;
00377       CurrentLevel[CurrentRow[Len-1]] = 0;
00378 
00379       // Merge List with rows in U
00380 
00381       First = CurrentRow[0];
00382       Next = First;
00383       while (Next < i) {
00384         int PrevInList = Next;
00385         int NextInList = LinkList[Next];
00386         int RowU = Next;
00387         // Get Indices for this row of U
00388         ArrayView<const local_ordinal_type> IndicesU;
00389         U_Graph_->getLocalRowView (RowU, IndicesU);
00390         // FIXME (mfh 23 Dec 2013) size() returns ptrdiff_t, not int.
00391         int LengthRowU = IndicesU.size ();
00392 
00393         int ii;
00394 
00395         // Scan RowU
00396 
00397         for (ii = 0; ii < LengthRowU; /*nop*/) {
00398           int CurInList = IndicesU[ii];
00399           if (CurInList < NextInList) {
00400             // new fill-in
00401             int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00402             if (NewLevel <= LevelFill_) {
00403               LinkList[PrevInList]  = CurInList;
00404               LinkList[CurInList] = NextInList;
00405               PrevInList = CurInList;
00406               CurrentLevel[CurInList] = NewLevel;
00407             }
00408             ii++;
00409           }
00410           else if (CurInList == NextInList) {
00411             PrevInList = NextInList;
00412             NextInList = LinkList[PrevInList];
00413             int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00414             CurrentLevel[CurInList] = std::min (CurrentLevel[CurInList], NewLevel);
00415             ii++;
00416           }
00417           else { // (CurInList > NextInList)
00418             PrevInList = NextInList;
00419             NextInList = LinkList[PrevInList];
00420           }
00421         }
00422         Next = LinkList[Next];
00423       }
00424 
00425       // Put pattern into L and U
00426 
00427       CurrentRow.resize (0);
00428 
00429       Next = First;
00430 
00431       // Lower
00432 
00433       while (Next < i) {
00434         CurrentRow.push_back (Next);
00435         Next = LinkList[Next];
00436       }
00437 
00438       // FIXME (mfh 23 Dec 2013) It's not clear to me that
00439       // removeLocalIndices works like people expect it to work.  In
00440       // particular, it does not actually change the column Map.
00441       L_Graph_->removeLocalIndices (i); // Delete current set of Indices
00442       if (CurrentRow.size() > 0) {
00443         L_Graph_->insertLocalIndices (i, CurrentRow ());
00444       }
00445 
00446       // Diagonal
00447 
00448       TEUCHOS_TEST_FOR_EXCEPTION(
00449         Next != i, std::runtime_error,
00450         "Ifpack2::IlukGraph::initialize: FATAL: U has zero diagonal")
00451 
00452       LevelsRowU[0] = CurrentLevel[Next];
00453       Next = LinkList[Next];
00454 
00455       // Upper
00456 
00457       CurrentRow.resize (0);
00458       LenU = 0;
00459 
00460       while (Next < NumMyRows) {
00461         LevelsRowU[LenU+1] = CurrentLevel[Next];
00462         CurrentRow.push_back (Next);
00463         ++LenU;
00464         Next = LinkList[Next];
00465       }
00466 
00467       // FIXME (mfh 23 Dec 2013) It's not clear to me that
00468       // removeLocalIndices works like people expect it to work.  In
00469       // particular, it does not actually change the column Map.
00470 
00471       U_Graph_->removeLocalIndices (i); // Delete current set of Indices
00472       if (LenU > 0) {
00473         U_Graph_->insertLocalIndices (i, CurrentRow ());
00474       }
00475 
00476       // Allocate and fill Level info for this row
00477       Levels[i] = std::vector<int> (LenU+1);
00478       for (size_t jj=0; jj<LenU+1; jj++) {
00479         Levels[i][jj] = LevelsRowU[jj];
00480       }
00481     }
00482   }
00483 
00484   // Complete Fill steps
00485   RCP<const map_type> L_DomainMap = OverlapGraph_->getRowMap ();
00486   RCP<const map_type> L_RangeMap  = Graph_->getRangeMap ();
00487   RCP<const map_type> U_DomainMap = Graph_->getDomainMap ();
00488   RCP<const map_type> U_RangeMap  = OverlapGraph_->getRowMap ();
00489   L_Graph_->fillComplete (L_DomainMap, L_RangeMap);//DoOptimizeStorage is default here...
00490   U_Graph_->fillComplete (U_DomainMap, U_RangeMap);//DoOptimizeStorage is default here...
00491 
00492   reduceAll<int, size_t> (* (L_DomainMap->getComm ()), REDUCE_SUM, 1,
00493                           &NumMyDiagonals_, &NumGlobalDiagonals_);
00494 }
00495 
00496 } // namespace Ifpack2
00497 
00498 #endif /* IFPACK2_ILUK_GRAPH_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends