Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
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
56namespace Tpetra {
57namespace Details {
58
59namespace Impl {
65 template<class LocalGraphType, class LocalMapType>
67 public:
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
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:
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_device_type = typename crs_graph_type::local_graph_device_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.getLocalGraphDevice (), rowMap->getLocalMap (), colMap->getLocalMap ());
132 Kokkos::parallel_reduce (policy_type (0, G.getLocalNumRows ()), 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
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 {
169 (! G.supportsRowViews (), std::logic_error, "Not implemented!");
170
171 typename ::Tpetra::RowGraph<LO, GO, NT>::local_inds_host_view_type
173 const LO lclNumRows = static_cast<LO> (G.getLocalNumRows ());
174
175 LO diagCount = 0;
176 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
177 G.getLocalRowView (lclRow, lclColInds);
178 const LO numEnt = static_cast<LO> (lclColInds.size ());
179 if (numEnt != 0) {
180 const LO lclDiagCol = colMap->getLocalElement (rowMap->getGlobalElement (lclRow));
181 // If it's not in the column Map, then there's no diagonal entry.
182 if (lclDiagCol != LOT::invalid ()) {
183 // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize
184 // for the sorted case.
185 for (LO k = 0; k < numEnt; ++k) {
186 if (lclDiagCol == lclColInds[k]) {
187 ++diagCount;
188 break; // don't count duplicate entries
189 }
190 } // for each columm index in lclRow
191 } // if lclDiagCol is valid
192 } // numEnt != 0
193 } // for each lclRow
194
195 return diagCount;
196 } // if-else
197 }
198
200 template<class LO, class GO, class NT>
201 typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
203 {
204 using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;
205
206 const auto rowMap = G.getRowMap ();
207 const auto colMap = G.getColMap ();
208 if (rowMap.get () == nullptr || colMap.get () == nullptr) {
209 return 0; // this process does not participate
210 }
211 else {
212 using inds_type = typename ::Tpetra::RowGraph<LO,GO,NT>::nonconst_local_inds_host_view_type;
213 inds_type lclColIndsBuf("lclColIndsBuf",G.getLocalMaxNumRowEntries());
214 const LO lclNumRows = static_cast<LO> (G.getLocalNumRows ());
215
216 LO diagCount = 0;
217 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
218 size_t numEntSizeT = G.getNumEntriesInLocalRow (lclRow);
219 const LO numEnt = static_cast<LO> (numEntSizeT);
220
221 inds_type lclColInds = Kokkos::subview(lclColIndsBuf,std::make_pair(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
249 {
250 const auto rowMap = G.getRowMap ();
251 if (rowMap.get () == nullptr) {
252 return 0; // this process does not participate
253 }
254 else {
255 typename ::Tpetra::RowGraph<LO,GO,NT>::global_inds_host_view_type
257 const LO lclNumRows = static_cast<LO> (G.getLocalNumRows ());
258
259 LO diagCount = 0;
260 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
261 const GO gblRow = rowMap->getGlobalElement (lclRow);
262 G.getGlobalRowView (gblRow, gblColInds);
263 const LO numEnt = static_cast<LO> (gblColInds.size ());
264 if (numEnt != 0) {
265 // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
266 // the sorted case.
267 for (LO k = 0; k < numEnt; ++k) {
268 if (gblRow == gblColInds[k]) {
269 ++diagCount;
270 break; // don't count duplicate entries
271 }
272 } // for each column index in lclRow
273 } // if numEnt != 0
274 } // for each lclRow
275
276 return diagCount;
277 } // if-else
278 }
279
281 template<class LO, class GO, class NT>
282 typename ::Tpetra::RowGraph<LO, GO, NT>::local_ordinal_type
284 {
285 using gids_type = typename ::Tpetra::RowGraph<LO,GO,NT>::nonconst_global_inds_host_view_type ;
286 const auto rowMap = G.getRowMap ();
287 if (rowMap.get () == nullptr) {
288 return 0; // this process does not participate
289 }
290 else {
292 const LO lclNumRows = static_cast<LO> (G.getLocalNumRows ());
293
294 LO diagCount = 0;
295 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
296 size_t numEntSizeT = G.getNumEntriesInLocalRow (lclRow);
297 const LO numEnt = static_cast<LO> (numEntSizeT);
298 if (static_cast<LO> (gblColIndsBuf.size ()) < numEnt) {
299 Kokkos::resize(gblColIndsBuf,numEnt);
300 }
301
302 gids_type gblColInds = Kokkos::subview(gblColIndsBuf,std::make_pair((LO)0, numEnt));
303 const GO gblRow = rowMap->getGlobalElement (lclRow);
304 G.getGlobalRowCopy (gblRow, gblColInds, numEntSizeT);
305
306 if (numEnt != 0) {
307 // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
308 // the sorted case.
309 for (LO k = 0; k < numEnt; ++k) {
310 if (gblRow == gblColInds[k]) {
311 ++diagCount;
312 break; // don't count duplicate entries
313 }
314 } // for each column index in lclRow
315 } // if numEnt != 0
316 } // for each lclRow
317
318 return diagCount;
319 } // if-else
320 }
321
328 template<class MatrixType>
330 static typename MatrixType::local_ordinal_type
331 getLocalNumDiags (const MatrixType& A)
332 {
333 using LO = typename MatrixType::local_ordinal_type;
334 using GO = typename MatrixType::global_ordinal_type;
335 using NT = typename MatrixType::node_type;
336 using row_graph_type = ::Tpetra::RowGraph<LO, GO, NT>;
337
338 auto G = A.getGraph ();
339 if (G.get () == nullptr) {
340 return 0;
341 }
342 else {
344 }
345 }
346 };
347
349 template<class LO, class GO, class NT>
350 struct GetLocalNumDiags< ::Tpetra::RowGraph<LO, GO, NT> > {
351 static LO
352 getLocalNumDiags (const ::Tpetra::RowGraph<LO, GO, NT>& G)
353 {
354 using crs_graph_type = ::Tpetra::CrsGraph<LO, GO, NT>;
355
356 const crs_graph_type* G_crs = dynamic_cast<const crs_graph_type*> (&G);
357 if (G_crs != nullptr && G_crs->isFillComplete ()) {
358 return countLocalNumDiagsInFillCompleteGraph (*G_crs);
359 }
360 else {
361 if (G.isLocallyIndexed ()) {
362 if (G.supportsRowViews ()) {
364 }
365 else {
367 }
368 }
369 else if (G.isGloballyIndexed ()) {
370 if (G.supportsRowViews ()) {
372 }
373 else {
375 }
376 }
377 else { // G is empty
378 return 0;
379 }
380 }
381 }
382 };
383
385 template<class LO, class GO, class NT>
386 struct GetLocalNumDiags< ::Tpetra::CrsGraph<LO, GO, NT> > {
387 static LO
388 getLocalNumDiags (const ::Tpetra::CrsGraph<LO, GO, NT>& G)
389 {
390 using row_graph_type = ::Tpetra::RowGraph<LO, GO, NT>;
392 }
393 };
394} // namespace Impl
395
398template<class CrsGraphType>
399typename CrsGraphType::local_ordinal_type
404
407template<class CrsGraphType>
408typename CrsGraphType::global_ordinal_type
410{
411 using GO = typename CrsGraphType::global_ordinal_type;
412
413 const auto map = G.getRowMap ();
414 if (map.get () == nullptr) {
415 return GO (0); // this process should not participate
416 }
417 else {
418 const auto comm = map->getComm ();
419 if (comm.get () == nullptr) {
420 return GO (0); // this process should not participate
421 }
422 else {
423 const GO lclNumDiags = static_cast<GO> (getLocalNumDiags (G));
424
425 using Teuchos::REDUCE_SUM;
426 using Teuchos::reduceAll;
427 using Teuchos::outArg;
428
429 GO gblNumDiags {0};
431 return gblNumDiags;
432 }
433 }
434}
435
436} // namespace Details
437} // namespace Tpetra
438
439#endif // TPETRA_DETAILS_GETNUMDIAGS_HPP
440
Import KokkosSparse::OrdinalTraits, a traits class for "invalid" (flag) values of integer types,...
typename::Tpetra::RowGraph< LO, GO, NT >::local_ordinal_type countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithoutRowViews(const ::Tpetra::RowGraph< LO, GO, NT > &G)
Return local number of diagonal entries.
typename::Tpetra::RowGraph< LO, GO, NT >::local_ordinal_type countLocalNumDiagsInNonFillCompleteGloballyIndexedGraphWithRowViews(const ::Tpetra::RowGraph< LO, GO, NT > &G)
Return local number of diagonal entries.
MapType::local_ordinal_type getLocalDiagonalColumnIndex(const typename MapType::local_ordinal_type lclRow, const MapType &rowMap, const MapType &colMap)
Local columm index of diagonal entry.
typename::Tpetra::RowGraph< LO, GO, NT >::local_ordinal_type countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithoutRowViews(const ::Tpetra::RowGraph< LO, GO, NT > &G)
Return local number of diagonal entries.
typename::Tpetra::RowGraph< LO, GO, NT >::local_ordinal_type countLocalNumDiagsInNonFillCompleteLocallyIndexedGraphWithRowViews(const ::Tpetra::RowGraph< LO, GO, NT > &G)
Return local number of diagonal entries.
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
Struct that holds views of the contents of a CrsMatrix.
Kokkos::parallel_reduce functor for counting the local number of diagonal entries in a sparse graph.
KOKKOS_INLINE_FUNCTION void operator()(const typename LocalMapType::local_ordinal_type lclRow, result_type &diagCount) const
Reduction function: result is (diagonal count, error count).
An abstract interface for graphs accessed by rows.
Implementation details of Tpetra.
CrsGraphType::global_ordinal_type getGlobalNumDiags(const CrsGraphType &G)
Number of populated diagonal entries in the given sparse graph, over all processes in the graph'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.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Implementation of Tpetra::Details::getLocalNumDiags (see below).