IFPACK  Development
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();
430  Epetra_CrsGraph & L = (Epetra_CrsGraph &) A.L_Graph();
431  Epetra_CrsGraph & U = (Epetra_CrsGraph &) A.U_Graph();
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 }
virtual const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
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...
virtual const Epetra_BlockMap & RangeMap() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
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.
virtual ~Ifpack_IlukGraph()
Ifpack_IlukGraph Destructor.
virtual int ConstructFilledGraph()
Does the actual construction of the graph.
friend std::ostream & operator<<(std::ostream &os, const Ifpack_IlukGraph &A)
<< operator will work for Ifpack_IlukGraph.