Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_ConstructLevelFillGraph.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_CONSTRUCTLEVELFILLGRAPH_HPP
00031 #define IFPACK2_CONSTRUCTLEVELFILLGRAPH_HPP
00032 
00033 #include "Ifpack2_ConfigDefs.hpp"
00034 #include <Teuchos_Describable.hpp>
00035 #include "Tpetra_RowGraph.hpp" 
00036 #include "Tpetra_CrsGraph.hpp"
00037 #include "Teuchos_RefCountPtr.hpp"
00038 
00039 namespace Teuchos {
00040   class ParameterList;
00041 }
00042 
00043 namespace Ifpack2 {
00044   
00045   
00047   
00070   template<class LocalOrdinal, class GlobalOrdinal = LocalOrdinal>
00071   void ConstructLevelFillGraph(const RowGraph<LocalOrdinal,GlobalOrdinal> & userGraph, 
00072                  Teuchos::RCP<CrsGraph<LocalOrdinal,GlobalOrdinal>> graphL, Teuchos::RCP<CrsGraph<LocalOrdinal,GlobalOrdinal>> graphU) {
00073   
00074   LocalOrdinal i, j;
00075   LocalOrdinal numIn, numL, numU;
00076   LocalOrdinal numMyDiagonals = 0;
00077   bool diagFound;
00078   
00079   
00080   
00081   graphL = Teuchos::rcp(new CrsGraph(userGraph.getRowMap(), userGraph.getRowMap(),  0));
00082   graphU = Teuchos::rcp(new CrsGraph(userGraph.getRowMap(), userGraph.getRowMap(),  0));
00083   
00084   
00085   // Get Maximun Row length
00086   int MaxNumIndices = userGraph.getMaxNumIndices();
00087   
00088   Teuchos::ArrayView<LocalOrdinal> userGraphRow;
00089   vector<LocalOrdinal> L(MaxNumIndices);
00090   vector<LocalOrdinal> U(MaxNumIndices);
00091   LocalOrdinal * pL = &L[0];
00092   LocalOrdinal * pU = &U[0];
00093   
00094   LocalOrdinal numMyRows = userGraph.getNumMyRows();
00095   
00096   // First we copy the user's graph into L and U, regardless of fill level
00097   
00098   for (i=0; i< numMyRows; ++i) {
00099     
00100     
00101     userGraph.extractMyRowConstView(i, userGraphRow); // Get Indices
00102     
00103     
00104     // Split into L and U (we don't assume that indices are ordered).
00105     
00106     numL = 0; 
00107     numU = 0; 
00108     diagFound = false;
00109     numIn = userGraphRow.size();
00110     
00111     for (j=0; j< numIn; ++j) {
00112       int k = userGraphRow[j];
00113       
00114       if (k<numMyRows) { // Ignore column elements that are not in the square matrix
00115         if (k==i) {DiagFound = true;} 
00116         else if (k < i) {pL[numL++] = k;}
00117         else {pU[numU++] = k;}
00118       }
00119     }
00120     
00121     // Check in things for this row of L and U
00122     
00123     if (diagFound) numMyDiagonals++;
00124     if (NumL>0) {graphL->insertMyIndices(i, ArrayView(pL, numL));}
00125     if (NumU>0) {graphU->insertMyIndices(i, ArrayView(pU, numU));}
00126   }
00127   
00128   //  if (LevelFill_ > 0) {
00129   //    
00130   //    // Complete Fill steps
00131   //    graphL->fillComplete(userGraph.getRowMap(), userGraph.getRangeMap());
00132   //    graphU->fillComplete(userGraph.getDomainMap(), userGraph.getRowMap());
00133   //    
00134   //    // At this point L_Graph and U_Graph are filled with the pattern of input graph, 
00135   //    // sorted and have redundant indices (if any) removed.  Indices are zero based.
00136   //    // LevelFill is greater than zero, so continue...
00137   //    
00138   //    vector<vector<LocalOrdinal> > levels(numMyRows);
00139   //    vector<LocalOrdinal> linkList(numMyRows);
00140   //    vector<LocalOrdinal> currentLevel(numMyRows);
00141   //    vector<LocalOrdinal> currentRow(numMyRows);
00142   //    vector<LocalOrdinal> levelsRowU(numMyRows);
00143   //    
00144   //    for (i=0; i<numMyRows; ++i) {
00145   //      int First, Next;
00146   //      
00147   //      // copy column indices of row into workspace and sort them
00148   //      
00149   //      int LenL = L_Graph_->NumMyIndices(i);
00150   //      int LenU = U_Graph_->NumMyIndices(i);
00151   //      int Len = LenL + LenU + 1;
00152   //      
00153   //      EPETRA_CHK_ERR(L_Graph_->ExtractMyRowCopy(i, LenL, LenL, &CurrentRow[0]));      // Get L Indices
00154   //      CurrentRow[LenL] = i;                                     // Put in Diagonal
00155   //      int ierr1 = 0;
00156   //      if (LenU) {
00157   //        // Get U Indices
00158   //        ierr1 = U_Graph_->ExtractMyRowCopy(i, LenU, LenU, &CurrentRow[LenL+1]);
00159   //      }
00160   //      if (ierr1!=0) {
00161   //        cout << "ierr1 = "<< ierr1 << endl;
00162   //        cout << "i = " << i << endl;
00163   //        cout << "NumMyBlockRows_ = " << graphU->NumMyBlockRows() << endl;
00164   //      }
00165   //        
00166   //      // Construct linked list for current row
00167   //        
00168   //      for (j=0; j<Len-1; j++) {
00169   //        LinkList[CurrentRow[j]] = CurrentRow[j+1];
00170   //        CurrentLevel[CurrentRow[j]] = 0;
00171   //      }
00172   //        
00173   //        LinkList[CurrentRow[Len-1]] = NumMyBlockRows_;
00174   //        CurrentLevel[CurrentRow[Len-1]] = 0;
00175   //        
00176   //        // Merge List with rows in U
00177   //        
00178   //        First = CurrentRow[0];
00179   //        Next = First;
00180   //        while (Next < i)
00181   //        {
00182   //          int PrevInList = Next;
00183   //          int NextInList = LinkList[Next];
00184   //          int RowU = Next;
00185   //          int LengthRowU;
00186   //          int * IndicesU;
00187   //          // Get Indices for this row of U
00188   //          EPETRA_CHK_ERR(U_Graph_->ExtractMyRowView(RowU, LengthRowU, IndicesU));
00189   //          
00190   //          int ii;
00191   //          
00192   //          // Scan RowU
00193   //          
00194   //          for (ii=0; ii<LengthRowU; /*nop*/) {
00195   //            int CurInList = IndicesU[ii];
00196   //            if (CurInList < NextInList) {
00197   //              // new fill-in
00198   //              int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00199   //              if (NewLevel <= LevelFill_)
00200   //              {
00201   //                LinkList[PrevInList]  = CurInList;
00202   //                LinkList[CurInList] = NextInList;
00203   //                PrevInList = CurInList;
00204   //                CurrentLevel[CurInList] = NewLevel;
00205   //              }
00206   //              ii++;
00207   //            }
00208   //            else if (CurInList == NextInList)
00209   //            {
00210   //              PrevInList = NextInList;
00211   //              NextInList = LinkList[PrevInList];
00212   //              int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00213   //              CurrentLevel[CurInList] = EPETRA_MIN(CurrentLevel[CurInList], NewLevel);
00214   //              ii++;
00215   //            }
00216   //            else // (CurInList > NextInList)
00217   //            {
00218   //              PrevInList = NextInList;
00219   //              NextInList = LinkList[PrevInList];
00220   //            }
00221   //          }
00222   //          Next = LinkList[Next];
00223   //        }
00224   //        
00225   //        // Put pattern into L and U
00226   //        
00227   //        LenL = 0;
00228   //        
00229   //        Next = First;
00230   //        
00231   //        // Lower
00232   //        
00233   //        while (Next < i) {    
00234   //          CurrentRow[LenL++] = Next;
00235   //          Next = LinkList[Next];
00236   //        }
00237   //        
00238   //        EPETRA_CHK_ERR(L_Graph_->RemoveMyIndices(i)); // Delete current set of Indices
00239   //        int ierr11 = L_Graph_->InsertMyIndices(i, LenL, &CurrentRow[0]);
00240   //        if (ierr11 < 0) EPETRA_CHK_ERR(ierr1);
00241   //        
00242   //        // Diagonal
00243   //        
00244   //        if (Next != i) return(-2); // Fatal:  U has zero diagonal.
00245   //        else {
00246   //          LevelsRowU[0] = CurrentLevel[Next];
00247   //          Next = LinkList[Next];
00248   //        }
00249   //        
00250   //        // Upper
00251   //        
00252   //        LenU = 0;
00253   //        
00254   //        while (Next < NumMyBlockRows_) // Should be "Next < NumMyBlockRows_"?
00255   //        {
00256   //          LevelsRowU[LenU+1] = CurrentLevel[Next];
00257   //          CurrentRow[LenU++] = Next;
00258   //          Next = LinkList[Next];
00259   //        }
00260   //        
00261   //        EPETRA_CHK_ERR(U_Graph_->RemoveMyIndices(i)); // Delete current set of Indices
00262   //        int ierr2 = U_Graph_->InsertMyIndices(i, LenU, &CurrentRow[0]);
00263   //        if (ierr2<0) EPETRA_CHK_ERR(ierr2);
00264   //        
00265   //        // Allocate and fill Level info for this row
00266   //        Levels[i] = vector<int>(LenU+1);
00267   //        for (int jj=0; jj<LenU+1; jj++) Levels[i][jj] = LevelsRowU[jj];
00268   //        
00269   //        }
00270   //        }    
00271   
00272   // Complete Fill steps
00273   graphL->fillComplete(userGraph.getRowMap(), userGraph.getRangeMap());
00274   graphU->fillComplete(userGraph.getDomainMap(), userGraph.getRowMap());
00275   
00276   
00277   // Optimize graph storage
00278   
00279   graphL->optimizeStorage();
00280   graphU->optimizeStorage();
00281   
00282   // Compute global quantities
00283   
00284   GlobalOrdinal numGlobalDiagonals = 0;
00285   
00286   graphL->getComm().sumAll(&numMyDiagonals, &numGlobalDiagonals, 1));
00287   
00288   return;
00289 }
00290 
00291 } // namespace Ifpack2
00292 
00293 #endif /* IFPACK2_CONSTRUCTLEVELFILLGRAPH_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends