49 #ifndef IFPACK2_ILUK_GRAPH_HPP 50 #define IFPACK2_ILUK_GRAPH_HPP 55 #include "KokkosSparse_spiluk.hpp" 57 #include <Ifpack2_ConfigDefs.hpp> 58 #include <Teuchos_ParameterList.hpp> 59 #include <Teuchos_CommHelpers.hpp> 60 #include <Tpetra_CrsGraph.hpp> 61 #include <Tpetra_Details_WrappedDualView.hpp> 62 #include <Tpetra_Import.hpp> 63 #include <Ifpack2_CreateOverlapGraph.hpp> 64 #include <Ifpack2_Parameters.hpp> 99 template<
class GraphType,
class KKHandleType>
102 typedef typename GraphType::local_ordinal_type local_ordinal_type;
103 typedef typename GraphType::global_ordinal_type global_ordinal_type;
104 typedef typename GraphType::node_type node_type;
107 typedef Tpetra::RowGraph<local_ordinal_type,
111 typedef Tpetra::CrsGraph<local_ordinal_type,
117 typedef typename crs_graph_type::nonconst_global_inds_host_view_type nonconst_global_inds_host_view_type;
118 typedef typename crs_graph_type::nonconst_local_inds_host_view_type nonconst_local_inds_host_view_type;
119 typedef typename crs_graph_type::global_inds_host_view_type global_inds_host_view_type;
120 typedef typename crs_graph_type::local_inds_host_view_type local_inds_host_view_type;
133 IlukGraph (
const Teuchos::RCP<const GraphType>& G,
135 const int levelOverlap,
136 const double overalloc = 2.);
146 void setParameters (
const Teuchos::ParameterList& parameterlist);
160 void initialize(
const Teuchos::RCP<KKHandleType>& KernelHandle);
180 return OverlapGraph_;
187 typedef typename GraphType::map_type map_type;
210 void constructOverlapGraph();
212 Teuchos::RCP<const GraphType> Graph_;
213 Teuchos::RCP<const crs_graph_type> OverlapGraph_;
216 const double Overalloc_;
217 Teuchos::RCP<crs_graph_type> L_Graph_;
218 Teuchos::RCP<crs_graph_type> U_Graph_;
219 size_t NumMyDiagonals_;
220 size_t NumGlobalDiagonals_;
224 template<
class GraphType,
class KKHandleType>
228 const int levelOverlap,
229 const double overalloc)
231 LevelFill_ (levelFill),
232 LevelOverlap_ (levelOverlap),
233 Overalloc_ (overalloc),
235 NumGlobalDiagonals_ (0)
237 TEUCHOS_TEST_FOR_EXCEPTION(Overalloc_ <= 1., std::runtime_error,
238 "Ifpack2::IlukGraph: FATAL: overalloc must be greater than 1.")
242 template<
class GraphType,
class KKHandleType>
247 template<
class GraphType,
class KKHandleType>
251 getParameter (parameterlist,
"iluk level-of-fill", LevelFill_);
252 getParameter (parameterlist,
"iluk level-of-overlap", LevelOverlap_);
256 template<
class GraphType,
class KKHandleType>
261 if (OverlapGraph_ == Teuchos::null) {
262 OverlapGraph_ = createOverlapGraph<GraphType> (Graph_, LevelOverlap_);
267 template<
class GraphType,
class KKHandleType>
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_->getNodeMaxNumRowEntries ();
287 const int NumMyRows = OverlapGraph_->getRowMap ()->getNodeNumElements ();
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::Experimental::ceil(static_cast<double>(RowMaxNumIndices)
309 * Kokkos::Experimental::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_);
564 template<
class GraphType,
class KKHandleType>
567 using Teuchos::Array;
568 using Teuchos::ArrayView;
571 using Teuchos::REDUCE_SUM;
572 using Teuchos::reduceAll;
574 typedef typename crs_graph_type::local_graph_device_type local_graph_device_type;
575 typedef typename local_graph_device_type::size_type size_type;
576 typedef typename local_graph_device_type::data_type data_type;
577 typedef typename local_graph_device_type::array_layout array_layout;
578 typedef typename local_graph_device_type::device_type device_type;
580 typedef typename Kokkos::View<size_type*, array_layout, device_type> lno_row_view_t;
581 typedef typename Kokkos::View<data_type*, array_layout, device_type> lno_nonzero_view_t;
583 constructOverlapGraph();
587 const int NumMyRows = OverlapGraph_->getRowMap()->getNodeNumElements();
588 auto localOverlapGraph = OverlapGraph_->getLocalGraphDevice();
590 if (KernelHandle->get_spiluk_handle()->get_nrows() <
static_cast<size_type
>(NumMyRows)) {
591 KernelHandle->get_spiluk_handle()->reset_handle(NumMyRows,
592 KernelHandle->get_spiluk_handle()->get_nnzL(),
593 KernelHandle->get_spiluk_handle()->get_nnzU());
596 lno_row_view_t L_row_map(
"L_row_map", NumMyRows + 1);
597 lno_nonzero_view_t L_entries(
"L_entries", KernelHandle->get_spiluk_handle()->get_nnzL());
598 lno_row_view_t U_row_map(
"U_row_map", NumMyRows + 1);
599 lno_nonzero_view_t U_entries(
"U_entries", KernelHandle->get_spiluk_handle()->get_nnzU());
603 symbolicError =
false;
605 KokkosSparse::Experimental::spiluk_symbolic( KernelHandle.getRawPtr(), LevelFill_,
606 localOverlapGraph.row_map, localOverlapGraph.entries,
607 L_row_map, L_entries, U_row_map, U_entries );
609 catch (std::runtime_error &e) {
610 symbolicError =
true;
611 data_type nnzL =
static_cast<data_type
>(Overalloc_)*L_entries.extent(0);
612 data_type nnzU =
static_cast<data_type
>(Overalloc_)*U_entries.extent(0);
613 KernelHandle->get_spiluk_handle()->reset_handle(NumMyRows, nnzL, nnzU);
614 Kokkos::resize(L_entries, KernelHandle->get_spiluk_handle()->get_nnzL());
615 Kokkos::resize(U_entries, KernelHandle->get_spiluk_handle()->get_nnzU());
617 const int localSymbolicError = symbolicError ? 1 : 0;
618 int globalSymbolicError = 0;
619 reduceAll (* (OverlapGraph_->getRowMap ()->getComm ()), REDUCE_SUM, 1,
620 &localSymbolicError, &globalSymbolicError);
621 symbolicError = globalSymbolicError > 0;
622 }
while (symbolicError);
624 Kokkos::resize(L_entries, KernelHandle->get_spiluk_handle()->get_nnzL());
625 Kokkos::resize(U_entries, KernelHandle->get_spiluk_handle()->get_nnzU());
627 RCP<Teuchos::ParameterList> params = Teuchos::parameterList ();
628 params->set (
"Optimize Storage",
false);
630 L_Graph_ = rcp (
new crs_graph_type (OverlapGraph_->getRowMap (),
631 OverlapGraph_->getRowMap (),
632 L_row_map, L_entries));
633 U_Graph_ = rcp (
new crs_graph_type (OverlapGraph_->getRowMap (),
634 OverlapGraph_->getRowMap (),
635 U_row_map, U_entries));
637 RCP<const map_type> L_DomainMap = OverlapGraph_->getRowMap ();
638 RCP<const map_type> L_RangeMap = Graph_->getRangeMap ();
639 RCP<const map_type> U_DomainMap = Graph_->getDomainMap ();
640 RCP<const map_type> U_RangeMap = OverlapGraph_->getRowMap ();
641 L_Graph_->fillComplete (L_DomainMap, L_RangeMap, params);
642 U_Graph_->fillComplete (U_DomainMap, U_RangeMap, params);
IlukGraph(const Teuchos::RCP< const GraphType > &G, const int levelFill, const int levelOverlap, const double overalloc=2.)
Constructor.
Definition: Ifpack2_IlukGraph.hpp:226
int getLevelFill() const
The level of fill used to construct this graph.
Definition: Ifpack2_IlukGraph.hpp:163
void getParameter(const Teuchos::ParameterList ¶ms, const std::string &name, T &value)
Set a value from a ParameterList if a parameter with the specified name exists.
Definition: Ifpack2_Parameters.hpp:59
virtual ~IlukGraph()
IlukGraph Destructor.
Definition: Ifpack2_IlukGraph.hpp:243
Tpetra::CrsGraph< local_ordinal_type, global_ordinal_type, node_type > crs_graph_type
Tpetra::CrsGraph specialization used by this class.
Definition: Ifpack2_IlukGraph.hpp:113
Teuchos::RCP< const crs_graph_type > getOverlapGraph() const
Returns the the overlapped graph.
Definition: Ifpack2_IlukGraph.hpp:179
Teuchos::RCP< crs_graph_type > getU_Graph() const
Returns the graph of upper triangle of the ILU(k) graph as a Tpetra::CrsGraph.
Definition: Ifpack2_IlukGraph.hpp:174
size_t getNumGlobalDiagonals() const
Returns the global number of diagonals in the ILU(k) graph.
Definition: Ifpack2_IlukGraph.hpp:184
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:100
int getLevelOverlap() const
The level of overlap used to construct this graph.
Definition: Ifpack2_IlukGraph.hpp:166
void initialize()
Set up the graph structure of the L and U factors.
Definition: Ifpack2_IlukGraph.hpp:268
Teuchos::RCP< crs_graph_type > getL_Graph() const
Returns the graph of lower triangle of the ILU(k) graph as a Tpetra::CrsGraph.
Definition: Ifpack2_IlukGraph.hpp:169
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
Tpetra::RowGraph< local_ordinal_type, global_ordinal_type, node_type > row_graph_type
Tpetra::RowGraph specialization used by this class.
Definition: Ifpack2_IlukGraph.hpp:109
void setParameters(const Teuchos::ParameterList ¶meterlist)
Set parameters.
Definition: Ifpack2_IlukGraph.hpp:249