IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_IlukGraph.cpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "Ifpack_IlukGraph.h"
44 #include "Epetra_Object.h"
45 #include "Epetra_Comm.h"
46 #include "Epetra_Import.h"
47 
48 #include <Teuchos_ParameterList.hpp>
49 #include <ifp_parameters.h>
50 
51 //==============================================================================
52 Ifpack_IlukGraph::Ifpack_IlukGraph(const Epetra_CrsGraph & Graph_in, int LevelFill_in, int LevelOverlap_in)
53  : Graph_(Graph_in),
54  DomainMap_(Graph_in.DomainMap()),
55  RangeMap_(Graph_in.RangeMap()),
56  Comm_(Graph_in.Comm()),
57  LevelFill_(LevelFill_in),
58  LevelOverlap_(LevelOverlap_in),
59  IndexBase_(Graph_in.IndexBase64()),
60  NumGlobalRows_(Graph_in.NumGlobalRows64()),
61  NumGlobalCols_(Graph_in.NumGlobalCols64()),
62  NumGlobalBlockRows_(Graph_in.NumGlobalBlockRows64()),
63  NumGlobalBlockCols_(Graph_in.NumGlobalBlockCols64()),
64  NumGlobalBlockDiagonals_(0),
65  NumGlobalNonzeros_(0),
66  NumGlobalEntries_(0),
67  NumMyBlockRows_(Graph_in.NumMyBlockRows()),
68  NumMyBlockCols_(Graph_in.NumMyBlockCols()),
69  NumMyRows_(Graph_in.NumMyRows()),
70  NumMyCols_(Graph_in.NumMyCols()),
71  NumMyBlockDiagonals_(0),
72  NumMyNonzeros_(0),
73  NumMyEntries_(0)
74 {
75 }
76 
77 //==============================================================================
79  : Graph_(Graph_in.Graph_),
80  DomainMap_(Graph_in.DomainMap()),
81  RangeMap_(Graph_in.RangeMap()),
82  Comm_(Graph_in.Comm()),
83  OverlapGraph_(Graph_in.OverlapGraph_),
84  OverlapRowMap_(Graph_in.OverlapRowMap_),
85  OverlapImporter_(Graph_in.OverlapImporter_),
86  LevelFill_(Graph_in.LevelFill_),
87  LevelOverlap_(Graph_in.LevelOverlap_),
88  IndexBase_(Graph_in.IndexBase_),
89  NumGlobalRows_(Graph_in.NumGlobalRows_),
90  NumGlobalCols_(Graph_in.NumGlobalCols_),
91  NumGlobalBlockRows_(Graph_in.NumGlobalBlockRows_),
92  NumGlobalBlockCols_(Graph_in.NumGlobalBlockCols_),
93  NumGlobalBlockDiagonals_(Graph_in.NumGlobalBlockDiagonals_),
94  NumGlobalNonzeros_(Graph_in.NumGlobalNonzeros_),
95  NumGlobalEntries_(Graph_in.NumGlobalEntries_),
96  NumMyBlockRows_(Graph_in.NumMyBlockRows_),
97  NumMyBlockCols_(Graph_in.NumMyBlockCols_),
98  NumMyRows_(Graph_in.NumMyRows_),
99  NumMyCols_(Graph_in.NumMyCols_),
100  NumMyBlockDiagonals_(Graph_in.NumMyBlockDiagonals_),
101  NumMyNonzeros_(Graph_in.NumMyNonzeros_),
102  NumMyEntries_(Graph_in.NumMyEntries_)
103 {
104  Epetra_CrsGraph & L_Graph_In = Graph_in.L_Graph();
105  Epetra_CrsGraph & U_Graph_In = Graph_in.U_Graph();
106  L_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(L_Graph_In) );
107  U_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(U_Graph_In) );
108 }
109 
110 //==============================================================================
112 {
113 }
114 
115 //==============================================================================
116 int Ifpack_IlukGraph::SetParameters(const Teuchos::ParameterList& parameterlist,
117  bool cerr_warning_if_unused)
118 {
119  Ifpack::param_struct params;
120  params.int_params[Ifpack::level_fill-FIRST_INT_PARAM] = LevelFill_;
121  params.int_params[Ifpack::level_overlap-FIRST_INT_PARAM] = LevelOverlap_;
122 
123  Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
124 
125  LevelFill_ = params.int_params[Ifpack::level_fill-FIRST_INT_PARAM];
126  LevelOverlap_ = params.int_params[Ifpack::level_overlap-FIRST_INT_PARAM];
127  return(0);
128 }
129 
130 //==============================================================================
132 
133  OverlapGraph_ = Teuchos::rcp( (Epetra_CrsGraph *) &Graph_, false );
134  OverlapRowMap_ = Teuchos::rcp( (Epetra_BlockMap *) &Graph_.RowMap(), false );
135 
136  if (LevelOverlap_==0 || !Graph_.DomainMap().DistributedGlobal()) return(0); // Nothing to do
137 
138  Teuchos::RefCountPtr<Epetra_CrsGraph> OldGraph;
139  Teuchos::RefCountPtr<Epetra_BlockMap> OldRowMap;
140  Epetra_BlockMap * DomainMap_tmp = (Epetra_BlockMap *) &Graph_.DomainMap();
141  Epetra_BlockMap * RangeMap_tmp = (Epetra_BlockMap *) &Graph_.RangeMap();
142  for (int level=1; level <= LevelOverlap_; level++) {
143  OldGraph = OverlapGraph_;
144  OldRowMap = OverlapRowMap_;
145 
146  OverlapImporter_ = Teuchos::rcp( (Epetra_Import *) OldGraph->Importer(), false );
147  OverlapRowMap_ = Teuchos::rcp( new Epetra_BlockMap(OverlapImporter_->TargetMap()) );
148 
149 
150  if (level<LevelOverlap_)
151  OverlapGraph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, *OverlapRowMap_, 0) );
152  else
153  // On last iteration, we want to filter out all columns except those that correspond
154  // to rows in the graph. This assures that our matrix is square
155  OverlapGraph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, *OverlapRowMap_, *OverlapRowMap_, 0) );
156 
157  EPETRA_CHK_ERR(OverlapGraph_->Import( Graph_, *OverlapImporter_, Insert));
158  if (level<LevelOverlap_) {
159  EPETRA_CHK_ERR(OverlapGraph_->FillComplete(*DomainMap_tmp, *RangeMap_tmp));
160  }
161  else {
162  // Copy last OverlapImporter because we will use it later
163  OverlapImporter_ = Teuchos::rcp( new Epetra_Import(*OverlapRowMap_, *DomainMap_tmp) );
164  EPETRA_CHK_ERR(OverlapGraph_->FillComplete(*DomainMap_tmp, *RangeMap_tmp));
165  }
166  }
167 
168  NumMyBlockRows_ = OverlapGraph_->NumMyBlockRows();
169  NumMyBlockCols_ = OverlapGraph_->NumMyBlockCols();
170  NumMyRows_ = OverlapGraph_->NumMyRows();
171  NumMyCols_ = OverlapGraph_->NumMyCols();
172 
173  return(0);
174 }
175 
176 //==============================================================================
178  using std::cout;
179  using std::endl;
180 
181  int ierr = 0;
182  int i, j;
183  int * In=0;
184  int NumIn, NumL, NumU;
185  bool DiagFound;
186 
187 
188  EPETRA_CHK_ERR(ConstructOverlapGraph());
189 
190  L_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, OverlapGraph_->RowMap(), OverlapGraph_->RowMap(), 0) );
191  U_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, OverlapGraph_->RowMap(), OverlapGraph_->RowMap(), 0));
192 
193 
194  // Get Maximun Row length
195  int MaxNumIndices = OverlapGraph_->MaxNumIndices();
196 
197  std::vector<int> L(MaxNumIndices);
198  std::vector<int> U(MaxNumIndices);
199 
200  // First we copy the user's graph into L and U, regardless of fill level
201 
202  for (i=0; i< NumMyBlockRows_; i++) {
203 
204 
205  OverlapGraph_->ExtractMyRowView(i, NumIn, In); // Get Indices
206 
207 
208  // Split into L and U (we don't assume that indices are ordered).
209 
210  NumL = 0;
211  NumU = 0;
212  DiagFound = false;
213 
214  for (j=0; j< NumIn; j++) {
215  int k = In[j];
216 
217  if (k<NumMyBlockRows_) { // Ignore column elements that are not in the square matrix
218 
219  if (k==i) DiagFound = true;
220 
221  else if (k < i) {
222  L[NumL] = k;
223  NumL++;
224  }
225  else {
226  U[NumU] = k;
227  NumU++;
228  }
229  }
230  }
231 
232  // Check in things for this row of L and U
233 
234  if (DiagFound) NumMyBlockDiagonals_++;
235  if (NumL) L_Graph_->InsertMyIndices(i, NumL, &L[0]);
236  if (NumU) U_Graph_->InsertMyIndices(i, NumU, &U[0]);
237 
238  }
239 
240  if (LevelFill_ > 0) {
241 
242  // Complete Fill steps
243  Epetra_BlockMap * L_DomainMap = (Epetra_BlockMap *) &OverlapGraph_->RowMap();
244  Epetra_BlockMap * L_RangeMap = (Epetra_BlockMap *) &Graph_.RangeMap();
245  Epetra_BlockMap * U_DomainMap = (Epetra_BlockMap *) &Graph_.DomainMap();
246  Epetra_BlockMap * U_RangeMap = (Epetra_BlockMap *) &OverlapGraph_->RowMap();
247  EPETRA_CHK_ERR(L_Graph_->FillComplete(*L_DomainMap, *L_RangeMap));
248  EPETRA_CHK_ERR(U_Graph_->FillComplete(*U_DomainMap, *U_RangeMap));
249 
250  // At this point L_Graph and U_Graph are filled with the pattern of input graph,
251  // sorted and have redundant indices (if any) removed. Indices are zero based.
252  // LevelFill is greater than zero, so continue...
253 
254  int MaxRC = NumMyBlockRows_;
255  std::vector<std::vector<int> > Levels(MaxRC);
256  std::vector<int> LinkList(MaxRC);
257  std::vector<int> CurrentLevel(MaxRC);
258  std::vector<int> CurrentRow(MaxRC);
259  std::vector<int> LevelsRowU(MaxRC);
260 
261  for (i=0; i<NumMyBlockRows_; i++)
262  {
263  int First, Next;
264 
265  // copy column indices of row into workspace and sort them
266 
267  int LenL = L_Graph_->NumMyIndices(i);
268  int LenU = U_Graph_->NumMyIndices(i);
269  int Len = LenL + LenU + 1;
270 
271  EPETRA_CHK_ERR(L_Graph_->ExtractMyRowCopy(i, LenL, LenL, &CurrentRow[0])); // Get L Indices
272  CurrentRow[LenL] = i; // Put in Diagonal
273  //EPETRA_CHK_ERR(U_Graph_->ExtractMyRowCopy(i, LenU, LenU, CurrentRow+LenL+1)); // Get U Indices
274  int ierr1 = 0;
275  if (LenU) {
276  // Get U Indices
277  ierr1 = U_Graph_->ExtractMyRowCopy(i, LenU, LenU, &CurrentRow[LenL+1]);
278  }
279  if (ierr1!=0) {
280  cout << "ierr1 = "<< ierr1 << endl;
281  cout << "i = " << i << endl;
282  cout << "NumMyBlockRows_ = " << U_Graph_->NumMyBlockRows() << endl;
283  }
284 
285  // Construct linked list for current row
286 
287  for (j=0; j<Len-1; j++) {
288  LinkList[CurrentRow[j]] = CurrentRow[j+1];
289  CurrentLevel[CurrentRow[j]] = 0;
290  }
291 
292  LinkList[CurrentRow[Len-1]] = NumMyBlockRows_;
293  CurrentLevel[CurrentRow[Len-1]] = 0;
294 
295  // Merge List with rows in U
296 
297  First = CurrentRow[0];
298  Next = First;
299  while (Next < i)
300  {
301  int PrevInList = Next;
302  int NextInList = LinkList[Next];
303  int RowU = Next;
304  int LengthRowU;
305  int * IndicesU;
306  // Get Indices for this row of U
307  EPETRA_CHK_ERR(U_Graph_->ExtractMyRowView(RowU, LengthRowU, IndicesU));
308 
309  int ii;
310 
311  // Scan RowU
312 
313  for (ii=0; ii<LengthRowU; /*nop*/)
314  {
315  int CurInList = IndicesU[ii];
316  if (CurInList < NextInList)
317  {
318  // new fill-in
319  int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
320  if (NewLevel <= LevelFill_)
321  {
322  LinkList[PrevInList] = CurInList;
323  LinkList[CurInList] = NextInList;
324  PrevInList = CurInList;
325  CurrentLevel[CurInList] = NewLevel;
326  }
327  ii++;
328  }
329  else if (CurInList == NextInList)
330  {
331  PrevInList = NextInList;
332  NextInList = LinkList[PrevInList];
333  int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
334  CurrentLevel[CurInList] = EPETRA_MIN(CurrentLevel[CurInList], NewLevel);
335  ii++;
336  }
337  else // (CurInList > NextInList)
338  {
339  PrevInList = NextInList;
340  NextInList = LinkList[PrevInList];
341  }
342  }
343  Next = LinkList[Next];
344  }
345 
346  // Put pattern into L and U
347 
348  LenL = 0;
349 
350  Next = First;
351 
352  // Lower
353 
354  while (Next < i) {
355  CurrentRow[LenL++] = Next;
356  Next = LinkList[Next];
357  }
358 
359  EPETRA_CHK_ERR(L_Graph_->RemoveMyIndices(i)); // Delete current set of Indices
360  int ierr11 = L_Graph_->InsertMyIndices(i, LenL, &CurrentRow[0]);
361  if (ierr11 < 0) EPETRA_CHK_ERR(ierr1);
362 
363  // Diagonal
364 
365  if (Next != i) return(-2); // Fatal: U has zero diagonal.
366  else {
367  LevelsRowU[0] = CurrentLevel[Next];
368  Next = LinkList[Next];
369  }
370 
371  // Upper
372 
373  LenU = 0;
374 
375  while (Next < NumMyBlockRows_) // Should be "Next < NumMyBlockRows_"?
376  {
377  LevelsRowU[LenU+1] = CurrentLevel[Next];
378  CurrentRow[LenU++] = Next;
379  Next = LinkList[Next];
380  }
381 
382  EPETRA_CHK_ERR(U_Graph_->RemoveMyIndices(i)); // Delete current set of Indices
383  int ierr2 = U_Graph_->InsertMyIndices(i, LenU, &CurrentRow[0]);
384  if (ierr2<0) EPETRA_CHK_ERR(ierr2);
385 
386  // Allocate and fill Level info for this row
387  Levels[i] = std::vector<int>(LenU+1);
388  for (int jj=0; jj<LenU+1; jj++) Levels[i][jj] = LevelsRowU[jj];
389 
390  }
391  }
392 
393  // Complete Fill steps
394  Epetra_BlockMap L_DomainMap = (Epetra_BlockMap) OverlapGraph_->RowMap();
395  Epetra_BlockMap L_RangeMap = (Epetra_BlockMap) Graph_.RangeMap();
396  Epetra_BlockMap U_DomainMap = (Epetra_BlockMap) Graph_.DomainMap();
397  Epetra_BlockMap U_RangeMap = (Epetra_BlockMap) OverlapGraph_->RowMap();
398  EPETRA_CHK_ERR(L_Graph_->FillComplete(L_DomainMap, L_RangeMap));
399  EPETRA_CHK_ERR(U_Graph_->FillComplete(U_DomainMap, U_RangeMap));
400 
401  // Optimize graph storage
402 
403  EPETRA_CHK_ERR(L_Graph_->OptimizeStorage());
404  EPETRA_CHK_ERR(U_Graph_->OptimizeStorage());
405 
406  // Compute global quantities
407 
408  NumGlobalBlockDiagonals_ = 0;
409  long long NumMyBlockDiagonals_LL = NumMyBlockDiagonals_;
410  EPETRA_CHK_ERR(L_Graph_->Comm().SumAll(&NumMyBlockDiagonals_LL, &NumGlobalBlockDiagonals_, 1));
411 
412  NumGlobalNonzeros_ = L_Graph_->NumGlobalNonzeros64()+U_Graph_->NumGlobalNonzeros64();
413  NumMyNonzeros_ = L_Graph_->NumMyNonzeros()+U_Graph_->NumMyNonzeros();
414  NumGlobalEntries_ = L_Graph_->NumGlobalEntries64()+U_Graph_->NumGlobalEntries64();
415  NumMyEntries_ = L_Graph_->NumMyEntries()+U_Graph_->NumMyEntries();
416  return(ierr);
417 }
418 //==========================================================================
419 
420 // Non-member functions
421 
422 std::ostream& operator << (std::ostream& os, const Ifpack_IlukGraph& A)
423 {
424  using std::endl;
425 
426 /* Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
427  Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
428  int oldp = os.precision(12); */
429  int LevelFill = A.LevelFill();
432  os.width(14);
433  os << " Level of Fill = "; os << LevelFill;
434  os << endl;
435 
436  os.width(14);
437  os << " Graph of L = ";
438  os << endl;
439  os << L; // Let Epetra_CrsGraph handle the rest.
440 
441  os.width(14);
442  os << " Graph of U = ";
443  os << endl;
444  os << U; // Let Epetra_CrsGraph handle the rest.
445 
446  // Reset os flags
447 
448 /* os.setf(olda,ios::adjustfield);
449  os.setf(oldf,ios::floatfield);
450  os.precision(oldp); */
451 
452  return os;
453 }
const Epetra_BlockMap & RangeMap() const
bool DistributedGlobal() const
Ifpack_IlukGraph(const Epetra_CrsGraph &Graph_in, int LevelFill_in, int LevelOverlap_in)
Ifpack_IlukGraph constuctor.
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
int SetParameters(const Teuchos::ParameterList &parameterlist, bool cerr_warning_if_unused=false)
Set parameters using Teuchos::ParameterList object.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
const Epetra_BlockMap & DomainMap() const
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
virtual int ConstructOverlapGraph()
Does the actual construction of the overlap matrix graph.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
const Epetra_BlockMap & RowMap() const
virtual ~Ifpack_IlukGraph()
Ifpack_IlukGraph Destructor.
virtual int ConstructFilledGraph()
Does the actual construction of the graph.