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