IFPACK Development
Ifpack_IlukGraph.cpp
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) 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 #include "Ifpack_IlukGraph.h"
00030 #include "Epetra_Object.h"
00031 #include "Epetra_Comm.h"
00032 #include "Epetra_Import.h"
00033 
00034 #include <Teuchos_ParameterList.hpp>
00035 #include <ifp_parameters.h>
00036 
00037 //==============================================================================
00038 Ifpack_IlukGraph::Ifpack_IlukGraph(const Epetra_CrsGraph & Graph_in, int LevelFill_in, int LevelOverlap_in)
00039   : Graph_(Graph_in),
00040     DomainMap_(Graph_in.DomainMap()),
00041     RangeMap_(Graph_in.RangeMap()),
00042     Comm_(Graph_in.Comm()),
00043     LevelFill_(LevelFill_in),
00044     LevelOverlap_(LevelOverlap_in),
00045     IndexBase_(Graph_in.IndexBase()),
00046     NumGlobalRows_(Graph_in.NumGlobalRows()),
00047     NumGlobalCols_(Graph_in.NumGlobalCols()),
00048     NumGlobalBlockRows_(Graph_in.NumGlobalBlockRows()),
00049     NumGlobalBlockCols_(Graph_in.NumGlobalBlockCols()),
00050     NumGlobalBlockDiagonals_(0),
00051     NumGlobalNonzeros_(0),
00052     NumGlobalEntries_(0),
00053     NumMyBlockRows_(Graph_in.NumMyBlockRows()),
00054     NumMyBlockCols_(Graph_in.NumMyBlockCols()),
00055     NumMyRows_(Graph_in.NumMyRows()),
00056     NumMyCols_(Graph_in.NumMyCols()),
00057     NumMyBlockDiagonals_(0),
00058     NumMyNonzeros_(0),
00059     NumMyEntries_(0)
00060 {
00061 }
00062 
00063 //==============================================================================
00064 Ifpack_IlukGraph::Ifpack_IlukGraph(const Ifpack_IlukGraph & Graph_in) 
00065   : Graph_(Graph_in.Graph_),
00066     DomainMap_(Graph_in.DomainMap()),
00067     RangeMap_(Graph_in.RangeMap()),
00068     Comm_(Graph_in.Comm()),
00069     OverlapGraph_(Graph_in.OverlapGraph_),
00070     OverlapRowMap_(Graph_in.OverlapRowMap_),
00071     OverlapImporter_(Graph_in.OverlapImporter_),
00072     LevelFill_(Graph_in.LevelFill_),
00073     LevelOverlap_(Graph_in.LevelOverlap_),
00074     IndexBase_(Graph_in.IndexBase_),
00075     NumGlobalRows_(Graph_in.NumGlobalRows_),
00076     NumGlobalCols_(Graph_in.NumGlobalCols_),
00077     NumGlobalBlockRows_(Graph_in.NumGlobalBlockRows_),
00078     NumGlobalBlockCols_(Graph_in.NumGlobalBlockCols_),
00079     NumGlobalBlockDiagonals_(Graph_in.NumGlobalBlockDiagonals_),
00080     NumGlobalNonzeros_(Graph_in.NumGlobalNonzeros_),
00081     NumGlobalEntries_(Graph_in.NumGlobalEntries_),
00082     NumMyBlockRows_(Graph_in.NumMyBlockRows_),
00083     NumMyBlockCols_(Graph_in.NumMyBlockCols_),
00084     NumMyRows_(Graph_in.NumMyRows_),
00085     NumMyCols_(Graph_in.NumMyCols_),
00086     NumMyBlockDiagonals_(Graph_in.NumMyBlockDiagonals_),
00087     NumMyNonzeros_(Graph_in.NumMyNonzeros_),
00088     NumMyEntries_(Graph_in.NumMyEntries_)
00089 {
00090   Epetra_CrsGraph & L_Graph_In = Graph_in.L_Graph();
00091   Epetra_CrsGraph & U_Graph_In = Graph_in.U_Graph();
00092   L_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(L_Graph_In) );
00093   U_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(U_Graph_In) );
00094 }
00095 
00096 //==============================================================================
00097 Ifpack_IlukGraph::~Ifpack_IlukGraph()
00098 {
00099 }
00100 
00101 //==============================================================================
00102 int Ifpack_IlukGraph::SetParameters(const Teuchos::ParameterList& parameterlist,
00103                                     bool cerr_warning_if_unused)
00104 {
00105   Ifpack::param_struct params;
00106   params.int_params[Ifpack::level_fill-FIRST_INT_PARAM] = LevelFill_;
00107   params.int_params[Ifpack::level_overlap-FIRST_INT_PARAM] = LevelOverlap_;
00108 
00109   Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
00110 
00111   LevelFill_ = params.int_params[Ifpack::level_fill-FIRST_INT_PARAM];
00112   LevelOverlap_ = params.int_params[Ifpack::level_overlap-FIRST_INT_PARAM];
00113   return(0);
00114 }
00115 
00116 //==============================================================================
00117 int Ifpack_IlukGraph::ConstructOverlapGraph() {
00118 
00119   OverlapGraph_ = Teuchos::rcp( (Epetra_CrsGraph *) &Graph_, false );
00120   OverlapRowMap_ = Teuchos::rcp( (Epetra_BlockMap *) &Graph_.RowMap(), false );
00121   
00122   if (LevelOverlap_==0 || !Graph_.DomainMap().DistributedGlobal()) return(0); // Nothing to do
00123 
00124   Teuchos::RefCountPtr<Epetra_CrsGraph> OldGraph;
00125   Teuchos::RefCountPtr<Epetra_BlockMap> OldRowMap;
00126   Epetra_BlockMap * DomainMap_tmp = (Epetra_BlockMap *) &Graph_.DomainMap();
00127   Epetra_BlockMap * RangeMap_tmp = (Epetra_BlockMap *) &Graph_.RangeMap();
00128   for (int level=1; level <= LevelOverlap_; level++) {
00129     OldGraph = OverlapGraph_; 
00130     OldRowMap = OverlapRowMap_;
00131 
00132     OverlapImporter_ = Teuchos::rcp( (Epetra_Import *) OldGraph->Importer(), false );
00133     OverlapRowMap_ = Teuchos::rcp( new Epetra_BlockMap(OverlapImporter_->TargetMap()) );
00134 
00135     
00136     if (level<LevelOverlap_)
00137       OverlapGraph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, *OverlapRowMap_, 0) );
00138     else
00139       // On last iteration, we want to filter out all columns except those that correspond
00140       // to rows in the graph.  This assures that our matrix is square
00141       OverlapGraph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, *OverlapRowMap_, *OverlapRowMap_, 0) );
00142 
00143     EPETRA_CHK_ERR(OverlapGraph_->Import( Graph_, *OverlapImporter_, Insert));
00144     if (level<LevelOverlap_) {
00145       EPETRA_CHK_ERR(OverlapGraph_->FillComplete(*DomainMap_tmp, *RangeMap_tmp));
00146     }
00147     else {
00148       // Copy last OverlapImporter because we will use it later
00149       OverlapImporter_ = Teuchos::rcp( new Epetra_Import(*OverlapRowMap_, *DomainMap_tmp) );
00150       EPETRA_CHK_ERR(OverlapGraph_->FillComplete(*DomainMap_tmp, *RangeMap_tmp));
00151     }
00152   }
00153 
00154     NumMyBlockRows_ = OverlapGraph_->NumMyBlockRows();
00155     NumMyBlockCols_ = OverlapGraph_->NumMyBlockCols();
00156     NumMyRows_ = OverlapGraph_->NumMyRows();
00157     NumMyCols_ = OverlapGraph_->NumMyCols();
00158 
00159   return(0);
00160 }
00161 
00162 //==============================================================================
00163 int Ifpack_IlukGraph::ConstructFilledGraph() {
00164   int ierr = 0;
00165   int i, j;
00166   int * In=0;
00167   int NumIn, NumL, NumU;
00168   bool DiagFound;
00169 
00170   
00171   EPETRA_CHK_ERR(ConstructOverlapGraph());
00172 
00173   L_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, OverlapGraph_->RowMap(), OverlapGraph_->RowMap(),  0) );
00174   U_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, OverlapGraph_->RowMap(), OverlapGraph_->RowMap(),  0));
00175 
00176 
00177   // Get Maximun Row length
00178   int MaxNumIndices = OverlapGraph_->MaxNumIndices();
00179 
00180   vector<int> L(MaxNumIndices);
00181   vector<int> U(MaxNumIndices);
00182 
00183   // First we copy the user's graph into L and U, regardless of fill level
00184 
00185   for (i=0; i< NumMyBlockRows_; i++) {
00186 
00187 
00188     OverlapGraph_->ExtractMyRowView(i, NumIn, In); // Get Indices
00189 
00190     
00191     // Split into L and U (we don't assume that indices are ordered).
00192     
00193     NumL = 0; 
00194     NumU = 0; 
00195     DiagFound = false;
00196     
00197     for (j=0; j< NumIn; j++) {
00198       int k = In[j];
00199 
00200       if (k<NumMyBlockRows_) { // Ignore column elements that are not in the square matrix
00201 
00202     if (k==i) DiagFound = true;
00203 
00204     else if (k < i) {
00205       L[NumL] = k;
00206       NumL++;
00207     }
00208     else {
00209       U[NumU] = k;
00210       NumU++;
00211     }
00212       }
00213     }
00214     
00215     // Check in things for this row of L and U
00216 
00217     if (DiagFound) NumMyBlockDiagonals_++;
00218     if (NumL) L_Graph_->InsertMyIndices(i, NumL, &L[0]);
00219     if (NumU) U_Graph_->InsertMyIndices(i, NumU, &U[0]);
00220     
00221   }
00222 
00223   if (LevelFill_ > 0) {
00224 
00225     // Complete Fill steps
00226     Epetra_BlockMap * L_DomainMap = (Epetra_BlockMap *) &OverlapGraph_->RowMap();
00227     Epetra_BlockMap * L_RangeMap = (Epetra_BlockMap *) &Graph_.RangeMap();
00228     Epetra_BlockMap * U_DomainMap = (Epetra_BlockMap *) &Graph_.DomainMap();
00229     Epetra_BlockMap * U_RangeMap = (Epetra_BlockMap *) &OverlapGraph_->RowMap();
00230     EPETRA_CHK_ERR(L_Graph_->FillComplete(*L_DomainMap, *L_RangeMap));
00231     EPETRA_CHK_ERR(U_Graph_->FillComplete(*U_DomainMap, *U_RangeMap));
00232 
00233     // At this point L_Graph and U_Graph are filled with the pattern of input graph, 
00234     // sorted and have redundant indices (if any) removed.  Indices are zero based.
00235     // LevelFill is greater than zero, so continue...
00236 
00237     int MaxRC = NumMyBlockRows_;
00238     vector<vector<int> > Levels(MaxRC);
00239     vector<int> LinkList(MaxRC);
00240     vector<int> CurrentLevel(MaxRC);
00241     vector<int> CurrentRow(MaxRC);
00242     vector<int> LevelsRowU(MaxRC);
00243 
00244     for (i=0; i<NumMyBlockRows_; i++)
00245     {
00246       int First, Next;
00247       
00248       // copy column indices of row into workspace and sort them
00249       
00250       int LenL = L_Graph_->NumMyIndices(i);
00251       int LenU = U_Graph_->NumMyIndices(i);
00252       int Len = LenL + LenU + 1;
00253       
00254       EPETRA_CHK_ERR(L_Graph_->ExtractMyRowCopy(i, LenL, LenL, &CurrentRow[0]));      // Get L Indices
00255       CurrentRow[LenL] = i;                                     // Put in Diagonal
00256       //EPETRA_CHK_ERR(U_Graph_->ExtractMyRowCopy(i, LenU, LenU, CurrentRow+LenL+1)); // Get U Indices
00257       int ierr1 = 0;
00258       if (LenU) {
00259         // Get U Indices
00260         ierr1 = U_Graph_->ExtractMyRowCopy(i, LenU, LenU, &CurrentRow[LenL+1]);
00261       }
00262       if (ierr1!=0) {
00263         cout << "ierr1 = "<< ierr1 << endl;
00264         cout << "i = " << i << endl;
00265         cout << "NumMyBlockRows_ = " << U_Graph_->NumMyBlockRows() << endl;
00266       }
00267       
00268       // Construct linked list for current row
00269       
00270       for (j=0; j<Len-1; j++) {
00271         LinkList[CurrentRow[j]] = CurrentRow[j+1];
00272         CurrentLevel[CurrentRow[j]] = 0;
00273       }
00274       
00275       LinkList[CurrentRow[Len-1]] = NumMyBlockRows_;
00276       CurrentLevel[CurrentRow[Len-1]] = 0;
00277       
00278       // Merge List with rows in U
00279       
00280       First = CurrentRow[0];
00281       Next = First;
00282       while (Next < i)
00283         {
00284       int PrevInList = Next;
00285       int NextInList = LinkList[Next];
00286       int RowU = Next;
00287       int LengthRowU;
00288       int * IndicesU;
00289       // Get Indices for this row of U
00290       EPETRA_CHK_ERR(U_Graph_->ExtractMyRowView(RowU, LengthRowU, IndicesU));
00291 
00292       int ii;
00293       
00294       // Scan RowU
00295       
00296       for (ii=0; ii<LengthRowU; /*nop*/)
00297             {
00298           int CurInList = IndicesU[ii];
00299           if (CurInList < NextInList)
00300                 {
00301           // new fill-in
00302           int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00303           if (NewLevel <= LevelFill_)
00304                     {
00305               LinkList[PrevInList]  = CurInList;
00306               LinkList[CurInList] = NextInList;
00307               PrevInList = CurInList;
00308               CurrentLevel[CurInList] = NewLevel;
00309                     }
00310           ii++;
00311                 }
00312           else if (CurInList == NextInList)
00313                 {
00314           PrevInList = NextInList;
00315           NextInList = LinkList[PrevInList];
00316           int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00317           CurrentLevel[CurInList] = EPETRA_MIN(CurrentLevel[CurInList], NewLevel);
00318           ii++;
00319                 }
00320           else // (CurInList > NextInList)
00321                 {
00322           PrevInList = NextInList;
00323           NextInList = LinkList[PrevInList];
00324                 }
00325             }
00326       Next = LinkList[Next];
00327         }
00328       
00329       // Put pattern into L and U
00330       
00331       LenL = 0;
00332 
00333       Next = First;
00334       
00335       // Lower
00336 
00337       while (Next < i) {      
00338     CurrentRow[LenL++] = Next;
00339     Next = LinkList[Next];
00340       }
00341 
00342       EPETRA_CHK_ERR(L_Graph_->RemoveMyIndices(i)); // Delete current set of Indices
00343       int ierr11 = L_Graph_->InsertMyIndices(i, LenL, &CurrentRow[0]);
00344       if (ierr11 < 0) EPETRA_CHK_ERR(ierr1);
00345 
00346       // Diagonal
00347 
00348       if (Next != i) return(-2); // Fatal:  U has zero diagonal.
00349       else {
00350     LevelsRowU[0] = CurrentLevel[Next];
00351     Next = LinkList[Next];
00352       }
00353 
00354       // Upper
00355 
00356       LenU = 0;
00357 
00358       while (Next < NumMyBlockRows_) // Should be "Next < NumMyBlockRows_"?
00359         {
00360       LevelsRowU[LenU+1] = CurrentLevel[Next];
00361       CurrentRow[LenU++] = Next;
00362       Next = LinkList[Next];
00363         }
00364 
00365       EPETRA_CHK_ERR(U_Graph_->RemoveMyIndices(i)); // Delete current set of Indices
00366       int ierr2 = U_Graph_->InsertMyIndices(i, LenU, &CurrentRow[0]);
00367       if (ierr2<0) EPETRA_CHK_ERR(ierr2);
00368 
00369       // Allocate and fill Level info for this row
00370       Levels[i] = vector<int>(LenU+1);
00371       for (int jj=0; jj<LenU+1; jj++) Levels[i][jj] = LevelsRowU[jj];
00372 
00373     }
00374   }    
00375 
00376   // Complete Fill steps
00377   Epetra_BlockMap L_DomainMap = (Epetra_BlockMap) OverlapGraph_->RowMap();
00378   Epetra_BlockMap L_RangeMap = (Epetra_BlockMap) Graph_.RangeMap();
00379   Epetra_BlockMap U_DomainMap = (Epetra_BlockMap) Graph_.DomainMap();
00380   Epetra_BlockMap U_RangeMap = (Epetra_BlockMap) OverlapGraph_->RowMap();
00381   EPETRA_CHK_ERR(L_Graph_->FillComplete(L_DomainMap, L_RangeMap));
00382   EPETRA_CHK_ERR(U_Graph_->FillComplete(U_DomainMap, U_RangeMap));
00383       
00384   // Optimize graph storage
00385   
00386   EPETRA_CHK_ERR(L_Graph_->OptimizeStorage());
00387   EPETRA_CHK_ERR(U_Graph_->OptimizeStorage());
00388 
00389   // Compute global quantities
00390 
00391   NumGlobalBlockDiagonals_ = 0;
00392 
00393   EPETRA_CHK_ERR(L_Graph_->Comm().SumAll(&NumMyBlockDiagonals_, &NumGlobalBlockDiagonals_, 1));
00394 
00395   NumGlobalNonzeros_ = L_Graph_->NumGlobalNonzeros()+U_Graph_->NumGlobalNonzeros();
00396   NumMyNonzeros_ = L_Graph_->NumMyNonzeros()+U_Graph_->NumMyNonzeros();
00397   NumGlobalEntries_ = L_Graph_->NumGlobalEntries()+U_Graph_->NumGlobalEntries();
00398   NumMyEntries_ = L_Graph_->NumMyEntries()+U_Graph_->NumMyEntries();
00399   return(ierr);
00400 }
00401 //==========================================================================
00402 
00403 // Non-member functions
00404 
00405 ostream& operator << (ostream& os, const Ifpack_IlukGraph& A)
00406 {
00407 /*  Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
00408   Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
00409   int oldp = os.precision(12); */
00410   int LevelFill = A.LevelFill();
00411   Epetra_CrsGraph & L = (Epetra_CrsGraph &) A.L_Graph();
00412   Epetra_CrsGraph & U = (Epetra_CrsGraph &) A.U_Graph();
00413   os.width(14);
00414   os <<  "     Level of Fill = "; os << LevelFill;
00415   os << endl;
00416 
00417   os.width(14);
00418   os <<  "     Graph of L = "; 
00419   os << endl;
00420   os << L; // Let Epetra_CrsGraph handle the rest.
00421 
00422   os.width(14);
00423   os <<  "     Graph of U = "; 
00424   os << endl;
00425   os << U; // Let Epetra_CrsGraph handle the rest.
00426  
00427   // Reset os flags
00428 
00429 /*  os.setf(olda,ios::adjustfield);
00430   os.setf(oldf,ios::floatfield);
00431   os.precision(oldp); */
00432 
00433   return os;
00434 }
 All Classes Files Functions Variables Enumerations Friends