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_)
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];
138 Teuchos::RefCountPtr<Epetra_CrsGraph> OldGraph;
139 Teuchos::RefCountPtr<Epetra_BlockMap> OldRowMap;
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_)
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) {
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];
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();
433 os <<
" Level of Fill = "; os << LevelFill;
437 os <<
" Graph of L = ";
442 os <<
" Graph of U = ";
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 ¶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...
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.