IFPACK Development
Loading...
Searching...
No Matches
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//==============================================================================
52Ifpack_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//==============================================================================
116int Ifpack_IlukGraph::SetParameters(const Teuchos::ParameterList& parameterlist,
117 bool cerr_warning_if_unused)
118{
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
422std::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}
Insert
Copy
bool DistributedGlobal() const
const Epetra_BlockMap & DomainMap() const
const Epetra_BlockMap & RowMap() const
const Epetra_BlockMap & RangeMap() const
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
int SetParameters(const Teuchos::ParameterList &parameterlist, bool cerr_warning_if_unused=false)
Set parameters using Teuchos::ParameterList object.
Ifpack_IlukGraph(const Epetra_CrsGraph &Graph_in, int LevelFill_in, int LevelOverlap_in)
Ifpack_IlukGraph constuctor.
virtual ~Ifpack_IlukGraph()
Ifpack_IlukGraph Destructor.
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 Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual int ConstructFilledGraph()
Does the actual construction of the graph.