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 = U_Graph_->ExtractMyRowCopy(i, LenU, LenU, &CurrentRow[LenL+1]); // Get U Indices
00258       if (ierr1!=0) {
00259     cout << "ierr1 = "<< ierr1 << endl;
00260     cout << "i = " << i << endl;
00261     cout << "NumMyBlockRows_ = " << U_Graph_->NumMyBlockRows() << endl;
00262       }
00263       
00264       // Construct linked list for current row
00265       
00266       for (j=0; j<Len-1; j++) {
00267     LinkList[CurrentRow[j]] = CurrentRow[j+1];
00268     CurrentLevel[CurrentRow[j]] = 0;
00269       }
00270       
00271       LinkList[CurrentRow[Len-1]] = NumMyBlockRows_;
00272       CurrentLevel[CurrentRow[Len-1]] = 0;
00273       
00274       // Merge List with rows in U
00275       
00276       First = CurrentRow[0];
00277       Next = First;
00278       while (Next < i)
00279         {
00280       int PrevInList = Next;
00281       int NextInList = LinkList[Next];
00282       int RowU = Next;
00283       int LengthRowU;
00284       int * IndicesU;
00285       // Get Indices for this row of U
00286       EPETRA_CHK_ERR(U_Graph_->ExtractMyRowView(RowU, LengthRowU, IndicesU));
00287 
00288       int ii;
00289       
00290       // Scan RowU
00291       
00292       for (ii=0; ii<LengthRowU; /*nop*/)
00293             {
00294           int CurInList = IndicesU[ii];
00295           if (CurInList < NextInList)
00296                 {
00297           // new fill-in
00298           int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00299           if (NewLevel <= LevelFill_)
00300                     {
00301               LinkList[PrevInList]  = CurInList;
00302               LinkList[CurInList] = NextInList;
00303               PrevInList = CurInList;
00304               CurrentLevel[CurInList] = NewLevel;
00305                     }
00306           ii++;
00307                 }
00308           else if (CurInList == NextInList)
00309                 {
00310           PrevInList = NextInList;
00311           NextInList = LinkList[PrevInList];
00312           int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
00313           CurrentLevel[CurInList] = EPETRA_MIN(CurrentLevel[CurInList], NewLevel);
00314           ii++;
00315                 }
00316           else // (CurInList > NextInList)
00317                 {
00318           PrevInList = NextInList;
00319           NextInList = LinkList[PrevInList];
00320                 }
00321             }
00322       Next = LinkList[Next];
00323         }
00324       
00325       // Put pattern into L and U
00326       
00327       LenL = 0;
00328 
00329       Next = First;
00330       
00331       // Lower
00332 
00333       while (Next < i) {      
00334     CurrentRow[LenL++] = Next;
00335     Next = LinkList[Next];
00336       }
00337 
00338       EPETRA_CHK_ERR(L_Graph_->RemoveMyIndices(i)); // Delete current set of Indices
00339       int ierr11 = L_Graph_->InsertMyIndices(i, LenL, &CurrentRow[0]);
00340       if (ierr11 < 0) EPETRA_CHK_ERR(ierr1);
00341 
00342       // Diagonal
00343 
00344       if (Next != i) return(-2); // Fatal:  U has zero diagonal.
00345       else {
00346     LevelsRowU[0] = CurrentLevel[Next];
00347     Next = LinkList[Next];
00348       }
00349        
00350 
00351       // Upper
00352 
00353       LenU = 0;
00354 
00355       while (Next < NumMyBlockRows_) // Should be "Next < NumMyBlockRows_"?
00356         {
00357       LevelsRowU[LenU+1] = CurrentLevel[Next];
00358       CurrentRow[LenU++] = Next;
00359       Next = LinkList[Next];
00360         }
00361 
00362       EPETRA_CHK_ERR(U_Graph_->RemoveMyIndices(i)); // Delete current set of Indices
00363       int ierr2 = U_Graph_->InsertMyIndices(i, LenU, &CurrentRow[0]);
00364       if (ierr2<0) EPETRA_CHK_ERR(ierr2);
00365 
00366       // Allocate and fill Level info for this row
00367       Levels[i] = vector<int>(LenU+1);
00368       for (int jj=0; jj<LenU+1; jj++) Levels[i][jj] = LevelsRowU[jj];
00369 
00370     }
00371   }    
00372 
00373   // Complete Fill steps
00374   Epetra_BlockMap L_DomainMap = (Epetra_BlockMap) OverlapGraph_->RowMap();
00375   Epetra_BlockMap L_RangeMap = (Epetra_BlockMap) Graph_.RangeMap();
00376   Epetra_BlockMap U_DomainMap = (Epetra_BlockMap) Graph_.DomainMap();
00377   Epetra_BlockMap U_RangeMap = (Epetra_BlockMap) OverlapGraph_->RowMap();
00378   EPETRA_CHK_ERR(L_Graph_->FillComplete(L_DomainMap, L_RangeMap));
00379   EPETRA_CHK_ERR(U_Graph_->FillComplete(U_DomainMap, U_RangeMap));
00380       
00381   // Optimize graph storage
00382   
00383   EPETRA_CHK_ERR(L_Graph_->OptimizeStorage());
00384   EPETRA_CHK_ERR(U_Graph_->OptimizeStorage());
00385 
00386   // Compute global quantities
00387 
00388   NumGlobalBlockDiagonals_ = 0;
00389 
00390   EPETRA_CHK_ERR(L_Graph_->Comm().SumAll(&NumMyBlockDiagonals_, &NumGlobalBlockDiagonals_, 1));
00391 
00392   NumGlobalNonzeros_ = L_Graph_->NumGlobalNonzeros()+U_Graph_->NumGlobalNonzeros();
00393   NumMyNonzeros_ = L_Graph_->NumMyNonzeros()+U_Graph_->NumMyNonzeros();
00394   NumGlobalEntries_ = L_Graph_->NumGlobalEntries()+U_Graph_->NumGlobalEntries();
00395   NumMyEntries_ = L_Graph_->NumMyEntries()+U_Graph_->NumMyEntries();
00396   return(ierr);
00397 }
00398 //==========================================================================
00399 
00400 // Non-member functions
00401 
00402 ostream& operator << (ostream& os, const Ifpack_IlukGraph& A)
00403 {
00404 /*  Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
00405   Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
00406   int oldp = os.precision(12); */
00407   int LevelFill = A.LevelFill();
00408   Epetra_CrsGraph & L = (Epetra_CrsGraph &) A.L_Graph();
00409   Epetra_CrsGraph & U = (Epetra_CrsGraph &) A.U_Graph();
00410   os.width(14);
00411   os <<  "     Level of Fill = "; os << LevelFill;
00412   os << endl;
00413 
00414   os.width(14);
00415   os <<  "     Graph of L = "; 
00416   os << endl;
00417   os << L; // Let Epetra_CrsGraph handle the rest.
00418 
00419   os.width(14);
00420   os <<  "     Graph of U = "; 
00421   os << endl;
00422   os << U; // Let Epetra_CrsGraph handle the rest.
00423  
00424   // Reset os flags
00425 
00426 /*  os.setf(olda,ios::adjustfield);
00427   os.setf(oldf,ios::floatfield);
00428   os.precision(oldp); */
00429 
00430   return os;
00431 }

Generated on Wed May 12 21:46:03 2010 for IFPACK by  doxygen 1.4.7