270 using Teuchos::Array;
271 using Teuchos::ArrayView;
274 using Teuchos::REDUCE_SUM;
275 using Teuchos::reduceAll;
277 size_t NumIn, NumL, NumU;
280 constructOverlapGraph();
283 const int MaxNumIndices = OverlapGraph_->getLocalMaxNumRowEntries ();
287 const int NumMyRows = OverlapGraph_->getRowMap ()->getLocalNumElements ();
289 using device_type =
typename node_type::device_type;
290 using execution_space =
typename device_type::execution_space;
291 using dual_view_type = Kokkos::DualView<size_t*,device_type>;
292 dual_view_type numEntPerRow_dv(
"numEntPerRow",NumMyRows);
293 Tpetra::Details::WrappedDualView<dual_view_type> numEntPerRow(numEntPerRow_dv);
295 const auto overalloc = Overalloc_;
296 const auto levelfill = LevelFill_;
299 auto numEntPerRow_d = numEntPerRow.getDeviceView(Tpetra::Access::OverwriteAll);
300 auto localOverlapGraph = OverlapGraph_->getLocalGraphDevice();
301 Kokkos::parallel_for(
"CountOverlapGraphRowEntries",
302 Kokkos::RangePolicy<execution_space>(0, NumMyRows),
303 KOKKOS_LAMBDA(
const int i)
306 int RowMaxNumIndices = localOverlapGraph.rowConst(i).length;
307 numEntPerRow_d(i) = (levelfill == 0) ? RowMaxNumIndices
308 : Kokkos::ceil(
static_cast<double>(RowMaxNumIndices)
309 * Kokkos::pow(overalloc, levelfill));
317 Teuchos::ArrayView<const size_t> a_numEntPerRow(numEntPerRow.getHostView(Tpetra::Access::ReadOnly).data(),NumMyRows);
319 OverlapGraph_->getRowMap (),
322 OverlapGraph_->getRowMap (),
325 Array<local_ordinal_type> L (MaxNumIndices);
326 Array<local_ordinal_type> U (MaxNumIndices);
332 for (
int i = 0; i< NumMyRows; ++i) {
333 local_inds_host_view_type my_indices;
334 OverlapGraph_->getLocalRowView (i, my_indices);
341 NumIn = my_indices.size();
343 for (
size_t j = 0; j < NumIn; ++j) {
344 const local_ordinal_type k = my_indices[j];
368 L_Graph_->insertLocalIndices (i, NumL, L.data());
371 U_Graph_->insertLocalIndices (i, NumU, U.data());
375 if (LevelFill_ > 0) {
377 RCP<const map_type> L_DomainMap = OverlapGraph_->getRowMap ();
378 RCP<const map_type> L_RangeMap = Graph_->getRangeMap ();
379 RCP<const map_type> U_DomainMap = Graph_->getDomainMap ();
380 RCP<const map_type> U_RangeMap = OverlapGraph_->getRowMap ();
381 RCP<Teuchos::ParameterList> params = Teuchos::parameterList ();
382 params->set (
"Optimize Storage",
false);
383 L_Graph_->fillComplete (L_DomainMap, L_RangeMap, params);
384 U_Graph_->fillComplete (U_DomainMap, U_RangeMap, params);
385 L_Graph_->resumeFill ();
386 U_Graph_->resumeFill ();
392 int MaxRC = NumMyRows;
393 std::vector<std::vector<int> > Levels(MaxRC);
394 std::vector<int> LinkList(MaxRC);
395 std::vector<int> CurrentLevel(MaxRC);
396 Array<local_ordinal_type> CurrentRow (MaxRC + 1);
397 std::vector<int> LevelsRowU(MaxRC);
400 for (
int i = 0; i < NumMyRows; ++i) {
405 size_t LenL = L_Graph_->getNumEntriesInLocalRow(i);
406 size_t LenU = U_Graph_->getNumEntriesInLocalRow(i);
407 size_t Len = LenL + LenU + 1;
408 CurrentRow.resize(Len);
409 nonconst_local_inds_host_view_type CurrentRow_view(CurrentRow.data(),CurrentRow.size());
410 L_Graph_->getLocalRowCopy(i, CurrentRow_view, LenL);
411 CurrentRow[LenL] = i;
413 ArrayView<local_ordinal_type> URowView = CurrentRow.view (LenL+1,LenU);
414 nonconst_local_inds_host_view_type URowView_v(URowView.data(),URowView.size());
417 U_Graph_->getLocalRowCopy (i, URowView_v, LenU);
422 for (
size_t j=0; j<Len-1; j++) {
423 LinkList[CurrentRow[j]] = CurrentRow[j+1];
424 CurrentLevel[CurrentRow[j]] = 0;
427 LinkList[CurrentRow[Len-1]] = NumMyRows;
428 CurrentLevel[CurrentRow[Len-1]] = 0;
432 First = CurrentRow[0];
435 int PrevInList = Next;
436 int NextInList = LinkList[Next];
439 local_inds_host_view_type IndicesU;
440 U_Graph_->getLocalRowView (RowU, IndicesU);
442 int LengthRowU = IndicesU.size ();
448 for (ii = 0; ii < LengthRowU; ) {
449 int CurInList = IndicesU[ii];
450 if (CurInList < NextInList) {
452 int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
453 if (NewLevel <= LevelFill_) {
454 LinkList[PrevInList] = CurInList;
455 LinkList[CurInList] = NextInList;
456 PrevInList = CurInList;
457 CurrentLevel[CurInList] = NewLevel;
461 else if (CurInList == NextInList) {
462 PrevInList = NextInList;
463 NextInList = LinkList[PrevInList];
464 int NewLevel = CurrentLevel[RowU] + Levels[RowU][ii+1] + 1;
465 CurrentLevel[CurInList] = std::min (CurrentLevel[CurInList],
470 PrevInList = NextInList;
471 NextInList = LinkList[PrevInList];
474 Next = LinkList[Next];
478 CurrentRow.resize(0);
484 CurrentRow.push_back(Next);
485 Next = LinkList[Next];
491 L_Graph_->removeLocalIndices (i);
492 if (CurrentRow.size() > 0) {
493 L_Graph_->insertLocalIndices (i, CurrentRow.size(),CurrentRow.data());
498 TEUCHOS_TEST_FOR_EXCEPTION(
499 Next != i, std::runtime_error,
500 "Ifpack2::IlukGraph::initialize: FATAL: U has zero diagonal")
502 LevelsRowU[0] = CurrentLevel[Next];
503 Next = LinkList[Next];
506 CurrentRow.resize(0);
509 while (Next < NumMyRows) {
510 LevelsRowU[LenU+1] = CurrentLevel[Next];
511 CurrentRow.push_back (Next);
513 Next = LinkList[Next];
520 U_Graph_->removeLocalIndices (i);
522 U_Graph_->insertLocalIndices (i, CurrentRow.size(),CurrentRow.data());
526 Levels[i] = std::vector<int> (LenU+1);
527 for (
size_t jj=0; jj<LenU+1; jj++) {
528 Levels[i][jj] = LevelsRowU[jj];
532 catch (std::runtime_error &e) {
534 auto numEntPerRow_d = numEntPerRow.getDeviceView(Tpetra::Access::OverwriteAll);
535 Kokkos::parallel_for(
"CountOverlapGraphRowEntries",
536 Kokkos::RangePolicy<execution_space>(0, NumMyRows),
537 KOKKOS_LAMBDA(
const int i)
539 const auto numRowEnt = numEntPerRow_d(i);
540 numEntPerRow_d(i) = ceil(
static_cast<double>((numRowEnt != 0 ? numRowEnt : 1)) * overalloc);
543 const int localInsertError = insertError ? 1 : 0;
544 int globalInsertError = 0;
545 reduceAll (* (OverlapGraph_->getRowMap ()->getComm ()), REDUCE_SUM, 1,
546 &localInsertError, &globalInsertError);
547 insertError = globalInsertError > 0;
549 }
while (insertError);
552 RCP<const map_type> L_DomainMap = OverlapGraph_->getRowMap ();
553 RCP<const map_type> L_RangeMap = Graph_->getRangeMap ();
554 RCP<const map_type> U_DomainMap = Graph_->getDomainMap ();
555 RCP<const map_type> U_RangeMap = OverlapGraph_->getRowMap ();
556 L_Graph_->fillComplete (L_DomainMap, L_RangeMap);
557 U_Graph_->fillComplete (U_DomainMap, U_RangeMap);
559 reduceAll<int, size_t> (* (L_DomainMap->getComm ()), REDUCE_SUM, 1,
560 &NumMyDiagonals_, &NumGlobalDiagonals_);