43 #include "Ifpack_IlukGraph.h" 44 #include "Epetra_Object.h" 45 #include "Epetra_Comm.h" 46 #include "Epetra_Import.h" 48 #include <Teuchos_ParameterList.hpp> 49 #include <ifp_parameters.h> 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),
67 NumMyBlockRows_(Graph_in.NumMyBlockRows()),
68 NumMyBlockCols_(Graph_in.NumMyBlockCols()),
69 NumMyRows_(Graph_in.NumMyRows()),
70 NumMyCols_(Graph_in.NumMyCols()),
71 NumMyBlockDiagonals_(0),
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_)
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) );
117 bool cerr_warning_if_unused)
120 params.int_params[Ifpack::level_fill-FIRST_INT_PARAM] = LevelFill_;
121 params.int_params[Ifpack::level_overlap-FIRST_INT_PARAM] = LevelOverlap_;
123 Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
125 LevelFill_ = params.int_params[Ifpack::level_fill-FIRST_INT_PARAM];
126 LevelOverlap_ = params.int_params[Ifpack::level_overlap-FIRST_INT_PARAM];
133 OverlapGraph_ = Teuchos::rcp( (Epetra_CrsGraph *) &Graph_,
false );
134 OverlapRowMap_ = Teuchos::rcp( (Epetra_BlockMap *) &Graph_.RowMap(), false );
136 if (LevelOverlap_==0 || !Graph_.DomainMap().DistributedGlobal())
return(0);
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_;
146 OverlapImporter_ = Teuchos::rcp( (Epetra_Import *) OldGraph->Importer(), false );
147 OverlapRowMap_ = Teuchos::rcp(
new Epetra_BlockMap(OverlapImporter_->TargetMap()) );
150 if (level<LevelOverlap_)
151 OverlapGraph_ = Teuchos::rcp(
new Epetra_CrsGraph(Copy, *OverlapRowMap_, 0) );
155 OverlapGraph_ = Teuchos::rcp(
new Epetra_CrsGraph(Copy, *OverlapRowMap_, *OverlapRowMap_, 0) );
157 EPETRA_CHK_ERR(OverlapGraph_->Import( Graph_, *OverlapImporter_, Insert));
158 if (level<LevelOverlap_) {
159 EPETRA_CHK_ERR(OverlapGraph_->FillComplete(*DomainMap_tmp, *RangeMap_tmp));
163 OverlapImporter_ = Teuchos::rcp(
new Epetra_Import(*OverlapRowMap_, *DomainMap_tmp) );
164 EPETRA_CHK_ERR(OverlapGraph_->FillComplete(*DomainMap_tmp, *RangeMap_tmp));
168 NumMyBlockRows_ = OverlapGraph_->NumMyBlockRows();
169 NumMyBlockCols_ = OverlapGraph_->NumMyBlockCols();
170 NumMyRows_ = OverlapGraph_->NumMyRows();
171 NumMyCols_ = OverlapGraph_->NumMyCols();
184 int NumIn, NumL, NumU;
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));
195 int MaxNumIndices = OverlapGraph_->MaxNumIndices();
197 std::vector<int> L(MaxNumIndices);
198 std::vector<int> U(MaxNumIndices);
202 for (i=0; i< NumMyBlockRows_; i++) {
205 OverlapGraph_->ExtractMyRowView(i, NumIn, In);
214 for (j=0; j< NumIn; j++) {
217 if (k<NumMyBlockRows_) {
219 if (k==i) DiagFound =
true;
234 if (DiagFound) NumMyBlockDiagonals_++;
235 if (NumL) L_Graph_->InsertMyIndices(i, NumL, &L[0]);
236 if (NumU) U_Graph_->InsertMyIndices(i, NumU, &U[0]);
240 if (LevelFill_ > 0) {
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));
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);
261 for (i=0; i<NumMyBlockRows_; i++)
267 int LenL = L_Graph_->NumMyIndices(i);
268 int LenU = U_Graph_->NumMyIndices(i);
269 int Len = LenL + LenU + 1;
271 EPETRA_CHK_ERR(L_Graph_->ExtractMyRowCopy(i, LenL, LenL, &CurrentRow[0]));
272 CurrentRow[LenL] = i;
277 ierr1 = U_Graph_->ExtractMyRowCopy(i, LenU, LenU, &CurrentRow[LenL+1]);
280 cout <<
"ierr1 = "<< ierr1 << endl;
281 cout <<
"i = " << i << endl;
282 cout <<
"NumMyBlockRows_ = " << U_Graph_->NumMyBlockRows() << endl;
287 for (j=0; j<Len-1; j++) {
288 LinkList[CurrentRow[j]] = CurrentRow[j+1];
289 CurrentLevel[CurrentRow[j]] = 0;
292 LinkList[CurrentRow[Len-1]] = NumMyBlockRows_;
293 CurrentLevel[CurrentRow[Len-1]] = 0;
297 First = CurrentRow[0];
301 int PrevInList = Next;
302 int NextInList = LinkList[Next];
307 EPETRA_CHK_ERR(U_Graph_->ExtractMyRowView(RowU, LengthRowU, IndicesU));
313 for (ii=0; ii<LengthRowU; )
315 int CurInList = IndicesU[ii];
316 if (CurInList < NextInList)
319 int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
320 if (NewLevel <= LevelFill_)
322 LinkList[PrevInList] = CurInList;
323 LinkList[CurInList] = NextInList;
324 PrevInList = CurInList;
325 CurrentLevel[CurInList] = NewLevel;
329 else if (CurInList == NextInList)
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);
339 PrevInList = NextInList;
340 NextInList = LinkList[PrevInList];
343 Next = LinkList[Next];
355 CurrentRow[LenL++] = Next;
356 Next = LinkList[Next];
359 EPETRA_CHK_ERR(L_Graph_->RemoveMyIndices(i));
360 int ierr11 = L_Graph_->InsertMyIndices(i, LenL, &CurrentRow[0]);
361 if (ierr11 < 0) EPETRA_CHK_ERR(ierr1);
365 if (Next != i)
return(-2);
367 LevelsRowU[0] = CurrentLevel[Next];
368 Next = LinkList[Next];
375 while (Next < NumMyBlockRows_)
377 LevelsRowU[LenU+1] = CurrentLevel[Next];
378 CurrentRow[LenU++] = Next;
379 Next = LinkList[Next];
382 EPETRA_CHK_ERR(U_Graph_->RemoveMyIndices(i));
383 int ierr2 = U_Graph_->InsertMyIndices(i, LenU, &CurrentRow[0]);
384 if (ierr2<0) EPETRA_CHK_ERR(ierr2);
387 Levels[i] = std::vector<int>(LenU+1);
388 for (
int jj=0; jj<LenU+1; jj++) Levels[i][jj] = LevelsRowU[jj];
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));
403 EPETRA_CHK_ERR(L_Graph_->OptimizeStorage());
404 EPETRA_CHK_ERR(U_Graph_->OptimizeStorage());
408 NumGlobalBlockDiagonals_ = 0;
409 long long NumMyBlockDiagonals_LL = NumMyBlockDiagonals_;
410 EPETRA_CHK_ERR(L_Graph_->Comm().SumAll(&NumMyBlockDiagonals_LL, &NumGlobalBlockDiagonals_, 1));
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();
430 Epetra_CrsGraph & L = (Epetra_CrsGraph &) A.
L_Graph();
431 Epetra_CrsGraph & U = (Epetra_CrsGraph &) A.
U_Graph();
433 os <<
" Level of Fill = "; os << LevelFill;
437 os <<
" Graph of L = ";
442 os <<
" Graph of U = ";
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 ¶meterlist, 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 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.