Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tpetra_Details_getNumDiags.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Tpetra: Templated Linear Algebra Services Package
4 //
5 // Copyright 2008 NTESS and the Tpetra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef TPETRA_DETAILS_GETNUMDIAGS_HPP
11 #define TPETRA_DETAILS_GETNUMDIAGS_HPP
12 
19 
20 #include "Tpetra_CrsGraph.hpp"
21 #include "Teuchos_CommHelpers.hpp"
23 
24 namespace Tpetra {
25 namespace Details {
26 
27 namespace Impl {
33 template <class LocalGraphType, class LocalMapType>
35  public:
36  CountLocalNumDiags(const LocalGraphType& G,
37  const LocalMapType& rowMap,
38  const LocalMapType& colMap)
39  : G_(G)
40  , rowMap_(rowMap)
41  , colMap_(colMap) {}
42 
43  // Result can't be more than the number of local rows, so
44  // local_ordinal_type is appropriate.
45  using result_type = typename LocalMapType::local_ordinal_type;
46 
48  KOKKOS_INLINE_FUNCTION void
49  operator()(const typename LocalMapType::local_ordinal_type lclRow,
50  result_type& diagCount) const {
51  using LO = typename LocalMapType::local_ordinal_type;
52  using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;
53 
54  auto G_row = G_.rowConst(lclRow);
55  const LO numEnt = G_row.length;
56  if (numEnt != 0) {
57  // Use global row and column indices to find the diagonal
58  // entry. Caller promises that local row index is in the row
59  // Map on the calling process.
60  const LO lclDiagCol = colMap_.getLocalElement(rowMap_.getGlobalElement(lclRow));
61  // If it's not in the column Map, then there's no diagonal entry.
62  if (lclDiagCol != LOT::invalid()) {
63  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
64  // the sorted case, but note that it requires operator[].
65  for (LO k = 0; k < numEnt; ++k) {
66  if (lclDiagCol == G_row(k)) {
67  ++diagCount;
68  break; // don't count duplicates
69  }
70  }
71  }
72  }
73  }
74 
75  private:
76  LocalGraphType G_;
77  LocalMapType rowMap_;
78  LocalMapType colMap_;
79 };
80 
81 template <class LO, class GO, class NT>
82 typename ::Tpetra::CrsGraph<LO, GO, NT>::local_ordinal_type
83 countLocalNumDiagsInFillCompleteGraph(const ::Tpetra::CrsGraph<LO, GO, NT>& G) {
84  using crs_graph_type = ::Tpetra::CrsGraph<LO, GO, NT>;
85  using local_map_type = typename crs_graph_type::map_type::local_map_type;
86  using local_graph_device_type = typename crs_graph_type::local_graph_device_type;
87  using functor_type = CountLocalNumDiags<local_graph_device_type, local_map_type>;
88  using execution_space = typename crs_graph_type::device_type::execution_space;
89  using policy_type = Kokkos::RangePolicy<execution_space, LO>;
90 
91  const auto rowMap = G.getRowMap();
92  const auto colMap = G.getColMap();
93  if (rowMap.get() == nullptr || colMap.get() == nullptr) {
94  return 0; // this process does not participate
95  } else {
96  LO lclNumDiags{0};
97  functor_type f(G.getLocalGraphDevice(), rowMap->getLocalMap(), colMap->getLocalMap());
98  Kokkos::parallel_reduce(policy_type(0, G.getLocalNumRows()), f, lclNumDiags);
99  return lclNumDiags;
100  }
101 }
102 
112 template <class MapType>
113 typename MapType::local_ordinal_type
114 getLocalDiagonalColumnIndex(const typename MapType::local_ordinal_type lclRow,
115  const MapType& rowMap,
116  const MapType& colMap) {
117  return colMap.getLocalElement(rowMap.getGlobalElement(lclRow));
118 }
119 
121 template <class LO, class GO, class NT>
122 typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
123 countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithRowViews(const ::Tpetra::RowGraph<LO, GO, NT>& G) {
124  using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;
125 
126  const auto rowMap = G.getRowMap();
127  const auto colMap = G.getColMap();
128  if (rowMap.get() == nullptr || colMap.get() == nullptr) {
129  return 0; // this process does not participate
130  } else {
131  TEUCHOS_TEST_FOR_EXCEPTION(!G.supportsRowViews(), std::logic_error, "Not implemented!");
132 
133  typename ::Tpetra::RowGraph<LO, GO, NT>::local_inds_host_view_type
134  lclColInds;
135  const LO lclNumRows = static_cast<LO>(G.getLocalNumRows());
136 
137  LO diagCount = 0;
138  for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
139  G.getLocalRowView(lclRow, lclColInds);
140  const LO numEnt = static_cast<LO>(lclColInds.size());
141  if (numEnt != 0) {
142  const LO lclDiagCol = colMap->getLocalElement(rowMap->getGlobalElement(lclRow));
143  // If it's not in the column Map, then there's no diagonal entry.
144  if (lclDiagCol != LOT::invalid()) {
145  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize
146  // for the sorted case.
147  for (LO k = 0; k < numEnt; ++k) {
148  if (lclDiagCol == lclColInds[k]) {
149  ++diagCount;
150  break; // don't count duplicate entries
151  }
152  } // for each columm index in lclRow
153  } // if lclDiagCol is valid
154  } // numEnt != 0
155  } // for each lclRow
156 
157  return diagCount;
158  } // if-else
159 }
160 
162 template <class LO, class GO, class NT>
163 typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
164 countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithoutRowViews(const ::Tpetra::RowGraph<LO, GO, NT>& G) {
165  using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;
166 
167  const auto rowMap = G.getRowMap();
168  const auto colMap = G.getColMap();
169  if (rowMap.get() == nullptr || colMap.get() == nullptr) {
170  return 0; // this process does not participate
171  } else {
172  using inds_type = typename ::Tpetra::RowGraph<LO, GO, NT>::nonconst_local_inds_host_view_type;
173  inds_type lclColIndsBuf("lclColIndsBuf", G.getLocalMaxNumRowEntries());
174  const LO lclNumRows = static_cast<LO>(G.getLocalNumRows());
175 
176  LO diagCount = 0;
177  for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
178  size_t numEntSizeT = G.getNumEntriesInLocalRow(lclRow);
179  const LO numEnt = static_cast<LO>(numEntSizeT);
180 
181  inds_type lclColInds = Kokkos::subview(lclColIndsBuf, std::make_pair(0, numEnt));
182  G.getLocalRowCopy(lclRow, lclColInds, numEntSizeT);
183 
184  if (numEnt != 0) {
185  const LO lclDiagCol =
186  colMap->getLocalElement(rowMap->getGlobalElement(lclRow));
187  // If it's not in the column Map, then there's no diagonal entry.
188  if (lclDiagCol != LOT::invalid()) {
189  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize
190  // for the sorted case.
191  for (LO k = 0; k < numEnt; ++k) {
192  if (lclDiagCol == lclColInds[k]) {
193  ++diagCount;
194  break; // don't count duplicate entries
195  }
196  } // for each columm index in lclRow
197  } // if lclDiagCol is valid
198  } // numEnt != 0
199  } // for each lclRow
200 
201  return diagCount;
202  } // if-else
203 }
204 
206 template <class LO, class GO, class NT>
207 typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
208 countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithRowViews(const ::Tpetra::RowGraph<LO, GO, NT>& G) {
209  const auto rowMap = G.getRowMap();
210  if (rowMap.get() == nullptr) {
211  return 0; // this process does not participate
212  } else {
213  typename ::Tpetra::RowGraph<LO, GO, NT>::global_inds_host_view_type
214  gblColInds;
215  const LO lclNumRows = static_cast<LO>(G.getLocalNumRows());
216 
217  LO diagCount = 0;
218  for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
219  const GO gblRow = rowMap->getGlobalElement(lclRow);
220  G.getGlobalRowView(gblRow, gblColInds);
221  const LO numEnt = static_cast<LO>(gblColInds.size());
222  if (numEnt != 0) {
223  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
224  // the sorted case.
225  for (LO k = 0; k < numEnt; ++k) {
226  if (gblRow == gblColInds[k]) {
227  ++diagCount;
228  break; // don't count duplicate entries
229  }
230  } // for each column index in lclRow
231  } // if numEnt != 0
232  } // for each lclRow
233 
234  return diagCount;
235  } // if-else
236 }
237 
239 template <class LO, class GO, class NT>
240 typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
241 countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithoutRowViews(const ::Tpetra::RowGraph<LO, GO, NT>& G) {
242  using gids_type = typename ::Tpetra::RowGraph<LO, GO, NT>::nonconst_global_inds_host_view_type;
243  const auto rowMap = G.getRowMap();
244  if (rowMap.get() == nullptr) {
245  return 0; // this process does not participate
246  } else {
247  gids_type gblColIndsBuf;
248  const LO lclNumRows = static_cast<LO>(G.getLocalNumRows());
249 
250  LO diagCount = 0;
251  for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
252  size_t numEntSizeT = G.getNumEntriesInLocalRow(lclRow);
253  const LO numEnt = static_cast<LO>(numEntSizeT);
254  if (static_cast<LO>(gblColIndsBuf.size()) < numEnt) {
255  Kokkos::resize(gblColIndsBuf, numEnt);
256  }
257 
258  gids_type gblColInds = Kokkos::subview(gblColIndsBuf, std::make_pair((LO)0, numEnt));
259  const GO gblRow = rowMap->getGlobalElement(lclRow);
260  G.getGlobalRowCopy(gblRow, gblColInds, numEntSizeT);
261 
262  if (numEnt != 0) {
263  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
264  // the sorted case.
265  for (LO k = 0; k < numEnt; ++k) {
266  if (gblRow == gblColInds[k]) {
267  ++diagCount;
268  break; // don't count duplicate entries
269  }
270  } // for each column index in lclRow
271  } // if numEnt != 0
272  } // for each lclRow
273 
274  return diagCount;
275  } // if-else
276 }
277 
284 template <class MatrixType>
286  static typename MatrixType::local_ordinal_type
287  getLocalNumDiags(const MatrixType& A) {
288  using LO = typename MatrixType::local_ordinal_type;
289  using GO = typename MatrixType::global_ordinal_type;
290  using NT = typename MatrixType::node_type;
291  using row_graph_type = ::Tpetra::RowGraph<LO, GO, NT>;
292 
293  auto G = A.getGraph();
294  if (G.get() == nullptr) {
295  return 0;
296  } else {
298  }
299  }
300 };
301 
303 template <class LO, class GO, class NT>
304 struct GetLocalNumDiags< ::Tpetra::RowGraph<LO, GO, NT> > {
305  static LO
306  getLocalNumDiags(const ::Tpetra::RowGraph<LO, GO, NT>& G) {
307  using crs_graph_type = ::Tpetra::CrsGraph<LO, GO, NT>;
308 
309  const crs_graph_type* G_crs = dynamic_cast<const crs_graph_type*>(&G);
310  if (G_crs != nullptr && G_crs->isFillComplete()) {
311  return countLocalNumDiagsInFillCompleteGraph(*G_crs);
312  } else {
313  if (G.isLocallyIndexed()) {
314  if (G.supportsRowViews()) {
315  return countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithRowViews(G);
316  } else {
317  return countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithoutRowViews(G);
318  }
319  } else if (G.isGloballyIndexed()) {
320  if (G.supportsRowViews()) {
321  return countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithRowViews(G);
322  } else {
323  return countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithoutRowViews(G);
324  }
325  } else { // G is empty
326  return 0;
327  }
328  }
329  }
330 };
331 
333 template <class LO, class GO, class NT>
334 struct GetLocalNumDiags< ::Tpetra::CrsGraph<LO, GO, NT> > {
335  static LO
336  getLocalNumDiags(const ::Tpetra::CrsGraph<LO, GO, NT>& G) {
337  using row_graph_type = ::Tpetra::RowGraph<LO, GO, NT>;
339  }
340 };
341 } // namespace Impl
342 
345 template <class CrsGraphType>
346 typename CrsGraphType::local_ordinal_type
347 getLocalNumDiags(const CrsGraphType& G) {
349 }
350 
353 template <class CrsGraphType>
354 typename CrsGraphType::global_ordinal_type
355 getGlobalNumDiags(const CrsGraphType& G) {
356  using GO = typename CrsGraphType::global_ordinal_type;
357 
358  const auto map = G.getRowMap();
359  if (map.get() == nullptr) {
360  return GO(0); // this process should not participate
361  } else {
362  const auto comm = map->getComm();
363  if (comm.get() == nullptr) {
364  return GO(0); // this process should not participate
365  } else {
366  const GO lclNumDiags = static_cast<GO>(getLocalNumDiags(G));
367 
368  using Teuchos::outArg;
369  using Teuchos::REDUCE_SUM;
370  using Teuchos::reduceAll;
371 
372  GO gblNumDiags{0};
373  reduceAll<int, GO>(*comm, REDUCE_SUM, lclNumDiags, outArg(gblNumDiags));
374  return gblNumDiags;
375  }
376  }
377 }
378 
379 } // namespace Details
380 } // namespace Tpetra
381 
382 #endif // TPETRA_DETAILS_GETNUMDIAGS_HPP
Teuchos::RCP< const map_type > getRowMap() const override
Returns the Map that describes the row distribution in this graph.
Import KokkosSparse::OrdinalTraits, a traits class for &quot;invalid&quot; (flag) values of integer types...
An abstract interface for graphs accessed by rows.
Kokkos::parallel_reduce functor for counting the local number of diagonal entries in a sparse graph...
CrsGraphType::global_ordinal_type getGlobalNumDiags(const CrsGraphType &G)
Number of populated diagonal entries in the given sparse graph, over all processes in the graph&#39;s (MP...
CrsGraphType::local_ordinal_type getLocalNumDiags(const CrsGraphType &G)
Number of populated diagonal entries in the given sparse graph, on the calling (MPI) process...
KOKKOS_INLINE_FUNCTION void operator()(const typename LocalMapType::local_ordinal_type lclRow, result_type &diagCount) const
Reduction function: result is (diagonal count, error count).
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
Implementation of Tpetra::Details::getLocalNumDiags (see below).