Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tpetra_Details_getNumDiags.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Tpetra: Templated Linear Algebra Services Package
5 // Copyright (2008) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 
42 #ifndef TPETRA_DETAILS_GETNUMDIAGS_HPP
43 #define TPETRA_DETAILS_GETNUMDIAGS_HPP
44 
51 
52 #include "Tpetra_CrsGraph.hpp"
53 #include "Teuchos_CommHelpers.hpp"
55 
56 namespace Tpetra {
57 namespace Details {
58 
59 namespace Impl {
65  template<class LocalGraphType, class LocalMapType>
67  public:
68  CountLocalNumDiags (const LocalGraphType& G,
69  const LocalMapType& rowMap,
70  const LocalMapType& colMap) :
71  G_ (G), rowMap_ (rowMap), colMap_ (colMap)
72  {}
73 
74  // Result can't be more than the number of local rows, so
75  // local_ordinal_type is appropriate.
76  using result_type = typename LocalMapType::local_ordinal_type;
77 
79  KOKKOS_INLINE_FUNCTION void
80  operator () (const typename LocalMapType::local_ordinal_type lclRow,
81  result_type& diagCount) const
82  {
83  using LO = typename LocalMapType::local_ordinal_type;
84  using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;
85 
86  auto G_row = G_.rowConst (lclRow);
87  const LO numEnt = G_row.length;
88  if (numEnt != 0) {
89  // Use global row and column indices to find the diagonal
90  // entry. Caller promises that local row index is in the row
91  // Map on the calling process.
92  const LO lclDiagCol = colMap_.getLocalElement (rowMap_.getGlobalElement (lclRow));
93  // If it's not in the column Map, then there's no diagonal entry.
94  if (lclDiagCol != LOT::invalid ()) {
95  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
96  // the sorted case, but note that it requires operator[].
97  for (LO k = 0; k < numEnt; ++k) {
98  if (lclDiagCol == G_row(k)) {
99  ++diagCount;
100  break; // don't count duplicates
101  }
102  }
103  }
104  }
105  }
106 
107  private:
108  LocalGraphType G_;
109  LocalMapType rowMap_;
110  LocalMapType colMap_;
111  };
112 
113  template<class LO, class GO, class NT>
114  typename ::Tpetra::CrsGraph<LO, GO, NT>::local_ordinal_type
115  countLocalNumDiagsInFillCompleteGraph (const ::Tpetra::CrsGraph<LO, GO, NT>& G)
116  {
117  using crs_graph_type = ::Tpetra::CrsGraph<LO, GO, NT>;
118  using local_map_type = typename crs_graph_type::map_type::local_map_type;
119  using local_graph_type = typename crs_graph_type::local_graph_type;
120  using functor_type = CountLocalNumDiags<local_graph_type, local_map_type>;
121  using execution_space = typename crs_graph_type::device_type::execution_space;
122  using policy_type = Kokkos::RangePolicy<execution_space, LO>;
123 
124  const auto rowMap = G.getRowMap ();
125  const auto colMap = G.getColMap ();
126  if (rowMap.get () == nullptr || colMap.get () == nullptr) {
127  return 0; // this process does not participate
128  }
129  else {
130  LO lclNumDiags {0};
131  functor_type f (G.getLocalGraph (), rowMap->getLocalMap (), colMap->getLocalMap ());
132  Kokkos::parallel_reduce (policy_type (0, G.getNodeNumRows ()), f, lclNumDiags);
133  return lclNumDiags;
134  }
135  }
136 
146  template<class MapType>
147  typename MapType::local_ordinal_type
148  getLocalDiagonalColumnIndex (const typename MapType::local_ordinal_type lclRow,
149  const MapType& rowMap,
150  const MapType& colMap)
151  {
152  return colMap.getLocalElement (rowMap.getGlobalElement (lclRow));
153  }
154 
156  template<class LO, class GO, class NT>
157  typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
158  countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithRowViews (const ::Tpetra::RowGraph<LO, GO, NT>& G)
159  {
160  using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;
161 
162  const auto rowMap = G.getRowMap ();
163  const auto colMap = G.getColMap ();
164  if (rowMap.get () == nullptr || colMap.get () == nullptr) {
165  return 0; // this process does not participate
166  }
167  else {
168  TEUCHOS_TEST_FOR_EXCEPTION
169  (! G.supportsRowViews (), std::logic_error, "Not implemented!");
170 
171  Teuchos::ArrayView<const LO> lclColInds;
172  const LO lclNumRows = static_cast<LO> (G.getNodeNumRows ());
173 
174  LO diagCount = 0;
175  for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
176  G.getLocalRowView (lclRow, lclColInds);
177  const LO numEnt = static_cast<LO> (lclColInds.size ());
178  if (numEnt != 0) {
179  const LO lclDiagCol = colMap->getLocalElement (rowMap->getGlobalElement (lclRow));
180  // If it's not in the column Map, then there's no diagonal entry.
181  if (lclDiagCol != LOT::invalid ()) {
182  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize
183  // for the sorted case.
184  for (LO k = 0; k < numEnt; ++k) {
185  if (lclDiagCol == lclColInds[k]) {
186  ++diagCount;
187  break; // don't count duplicate entries
188  }
189  } // for each columm index in lclRow
190  } // if lclDiagCol is valid
191  } // numEnt != 0
192  } // for each lclRow
193 
194  return diagCount;
195  } // if-else
196  }
197 
199  template<class LO, class GO, class NT>
200  typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
201  countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithoutRowViews (const ::Tpetra::RowGraph<LO, GO, NT>& G)
202  {
203  using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;
204 
205  const auto rowMap = G.getRowMap ();
206  const auto colMap = G.getColMap ();
207  if (rowMap.get () == nullptr || colMap.get () == nullptr) {
208  return 0; // this process does not participate
209  }
210  else {
211  Teuchos::Array<LO> lclColIndsBuf;
212  const LO lclNumRows = static_cast<LO> (G.getNodeNumRows ());
213 
214  LO diagCount = 0;
215  for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
216  size_t numEntSizeT = G.getNumEntriesInLocalRow (lclRow);
217  const LO numEnt = static_cast<LO> (numEntSizeT);
218  if (static_cast<LO> (lclColIndsBuf.size ()) < numEnt) {
219  lclColIndsBuf.resize (numEnt);
220  }
221  Teuchos::ArrayView<LO> lclColInds = lclColIndsBuf (0, numEnt);
222  G.getLocalRowCopy (lclRow, lclColInds, numEntSizeT);
223 
224  if (numEnt != 0) {
225  const LO lclDiagCol =
226  colMap->getLocalElement (rowMap->getGlobalElement (lclRow));
227  // If it's not in the column Map, then there's no diagonal entry.
228  if (lclDiagCol != LOT::invalid ()) {
229  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize
230  // for the sorted case.
231  for (LO k = 0; k < numEnt; ++k) {
232  if (lclDiagCol == lclColInds[k]) {
233  ++diagCount;
234  break; // don't count duplicate entries
235  }
236  } // for each columm index in lclRow
237  } // if lclDiagCol is valid
238  } // numEnt != 0
239  } // for each lclRow
240 
241  return diagCount;
242  } // if-else
243  }
244 
246  template<class LO, class GO, class NT>
247  typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
248  countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithRowViews (const ::Tpetra::RowGraph<LO, GO, NT>& G)
249  {
250  const auto rowMap = G.getRowMap ();
251  if (rowMap.get () == nullptr) {
252  return 0; // this process does not participate
253  }
254  else {
255  Teuchos::ArrayView<const GO> gblColInds;
256  const LO lclNumRows = static_cast<LO> (G.getNodeNumRows ());
257 
258  LO diagCount = 0;
259  for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
260  const GO gblRow = rowMap->getGlobalElement (lclRow);
261  G.getGlobalRowView (gblRow, gblColInds);
262  const LO numEnt = static_cast<LO> (gblColInds.size ());
263  if (numEnt != 0) {
264  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
265  // the sorted case.
266  for (LO k = 0; k < numEnt; ++k) {
267  if (gblRow == gblColInds[k]) {
268  ++diagCount;
269  break; // don't count duplicate entries
270  }
271  } // for each column index in lclRow
272  } // if numEnt != 0
273  } // for each lclRow
274 
275  return diagCount;
276  } // if-else
277  }
278 
280  template<class LO, class GO, class NT>
281  typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
282  countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithoutRowViews (const ::Tpetra::RowGraph<LO, GO, NT>& G)
283  {
284  const auto rowMap = G.getRowMap ();
285  if (rowMap.get () == nullptr) {
286  return 0; // this process does not participate
287  }
288  else {
289  Teuchos::Array<GO> gblColIndsBuf;
290  const LO lclNumRows = static_cast<LO> (G.getNodeNumRows ());
291 
292  LO diagCount = 0;
293  for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
294  size_t numEntSizeT = G.getNumEntriesInLocalRow (lclRow);
295  const LO numEnt = static_cast<LO> (numEntSizeT);
296  if (static_cast<LO> (gblColIndsBuf.size ()) < numEnt) {
297  gblColIndsBuf.resize (numEnt);
298  }
299  Teuchos::ArrayView<GO> gblColInds = gblColIndsBuf (0, numEnt);
300  const GO gblRow = rowMap->getGlobalElement (lclRow);
301  G.getGlobalRowCopy (gblRow, gblColInds, numEntSizeT);
302 
303  if (numEnt != 0) {
304  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
305  // the sorted case.
306  for (LO k = 0; k < numEnt; ++k) {
307  if (gblRow == gblColInds[k]) {
308  ++diagCount;
309  break; // don't count duplicate entries
310  }
311  } // for each column index in lclRow
312  } // if numEnt != 0
313  } // for each lclRow
314 
315  return diagCount;
316  } // if-else
317  }
318 
325  template<class MatrixType>
327  static typename MatrixType::local_ordinal_type
328  getLocalNumDiags (const MatrixType& A)
329  {
330  using LO = typename MatrixType::local_ordinal_type;
331  using GO = typename MatrixType::global_ordinal_type;
332  using NT = typename MatrixType::node_type;
333  using row_graph_type = ::Tpetra::RowGraph<LO, GO, NT>;
334 
335  auto G = A.getGraph ();
336  if (G.get () == nullptr) {
337  return 0;
338  }
339  else {
341  }
342  }
343  };
344 
346  template<class LO, class GO, class NT>
347  struct GetLocalNumDiags< ::Tpetra::RowGraph<LO, GO, NT> > {
348  static LO
349  getLocalNumDiags (const ::Tpetra::RowGraph<LO, GO, NT>& G)
350  {
351  using crs_graph_type = ::Tpetra::CrsGraph<LO, GO, NT>;
352 
353  const crs_graph_type* G_crs = dynamic_cast<const crs_graph_type*> (&G);
354  if (G_crs != nullptr && G_crs->isFillComplete ()) {
355  return countLocalNumDiagsInFillCompleteGraph (*G_crs);
356  }
357  else {
358  if (G.isLocallyIndexed ()) {
359  if (G.supportsRowViews ()) {
360  return countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithRowViews (G);
361  }
362  else {
363  return countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithoutRowViews (G);
364  }
365  }
366  else if (G.isGloballyIndexed ()) {
367  if (G.supportsRowViews ()) {
368  return countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithRowViews (G);
369  }
370  else {
371  return countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithoutRowViews (G);
372  }
373  }
374  else { // G is empty
375  return 0;
376  }
377  }
378  }
379  };
380 
382  template<class LO, class GO, class NT>
383  struct GetLocalNumDiags< ::Tpetra::CrsGraph<LO, GO, NT> > {
384  static LO
385  getLocalNumDiags (const ::Tpetra::CrsGraph<LO, GO, NT>& G)
386  {
387  using row_graph_type = ::Tpetra::RowGraph<LO, GO, NT>;
389  }
390  };
391 } // namespace Impl
392 
395 template<class CrsGraphType>
396 typename CrsGraphType::local_ordinal_type
397 getLocalNumDiags (const CrsGraphType& G)
398 {
400 }
401 
404 template<class CrsGraphType>
405 typename CrsGraphType::global_ordinal_type
406 getGlobalNumDiags (const CrsGraphType& G)
407 {
408  using GO = typename CrsGraphType::global_ordinal_type;
409 
410  const auto map = G.getRowMap ();
411  if (map.get () == nullptr) {
412  return GO (0); // this process should not participate
413  }
414  else {
415  const auto comm = map->getComm ();
416  if (comm.get () == nullptr) {
417  return GO (0); // this process should not participate
418  }
419  else {
420  const GO lclNumDiags = static_cast<GO> (getLocalNumDiags (G));
421 
422  using Teuchos::REDUCE_SUM;
423  using Teuchos::reduceAll;
424  using Teuchos::outArg;
425 
426  GO gblNumDiags {0};
427  reduceAll<int, GO> (*comm, REDUCE_SUM, lclNumDiags, outArg (gblNumDiags));
428  return gblNumDiags;
429  }
430  }
431 }
432 
433 } // namespace Details
434 } // namespace Tpetra
435 
436 #endif // TPETRA_DETAILS_GETNUMDIAGS_HPP
437 
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).