Ifpack_IlukGraph.cpp

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

Generated on Thu Sep 18 12:37:21 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1