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