Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tpetra_Details_determineLocalTriangularStructure.hpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Tpetra: Templated Linear Algebra Services Package
6 // Copyright (2008) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 */
43 
44 #ifndef TPETRA_DETAILS_DETERMINELOCALTRIANGULARSTRUCTURE_HPP
45 #define TPETRA_DETAILS_DETERMINELOCALTRIANGULARSTRUCTURE_HPP
46 
53 
54 #include "Kokkos_Core.hpp"
56 
57 namespace Tpetra {
58 namespace Details {
59 
63 template<class LO>
73 };
74 
75 namespace Impl {
89  template<class LocalGraphType, class LocalMapType>
91  public:
92  // Result can't be more than the number of local rows, so
93  // local_ordinal_type is appropriate.
94  using result_type =
96 
106  DetermineLocalTriangularStructure (const LocalGraphType& G,
107  const LocalMapType& rowMap,
108  const LocalMapType& colMap,
109  const bool ignoreMapsForTriangularStructure) :
110  G_ (G),
111  rowMap_ (rowMap),
112  colMap_ (colMap),
113  ignoreMapsForTriangularStructure_ (ignoreMapsForTriangularStructure)
114  {}
115 
117  KOKKOS_INLINE_FUNCTION void init (result_type& dst) const
118  {
119  dst.diagCount = 0;
120  dst.maxNumRowEnt = 0;
121  dst.couldBeLowerTriangular = true; // well, we don't know yet, do we?
122  dst.couldBeUpperTriangular = true; // ditto
123  }
124 
125  KOKKOS_INLINE_FUNCTION void
126  join (volatile result_type& dst,
127  const volatile result_type& src) const
128  {
129  dst.diagCount += src.diagCount;
130  dst.maxNumRowEnt = (src.maxNumRowEnt > dst.maxNumRowEnt) ?
131  src.maxNumRowEnt : dst.maxNumRowEnt;
132  dst.couldBeLowerTriangular &= src.couldBeLowerTriangular;
133  dst.couldBeUpperTriangular &= src.couldBeUpperTriangular;
134  }
135 
137  KOKKOS_INLINE_FUNCTION void
138  operator () (const typename LocalMapType::local_ordinal_type lclRow,
139  result_type& result) const
140  {
141  using LO = typename LocalMapType::local_ordinal_type;
142  using GO = typename LocalMapType::global_ordinal_type;
143  using LOT = typename ::Tpetra::Details::OrdinalTraits<LO>;
144 
145  auto G_row = G_.rowConst (lclRow);
146  const LO numEnt = G_row.length;
147  if (numEnt != 0) {
148  result.maxNumRowEnt = (numEnt > result.maxNumRowEnt) ?
149  numEnt : result.maxNumRowEnt;
150  // Use global row and column indices to find the diagonal
151  // entry. Caller promises that local row index is in the row
152  // Map on the calling process.
153  const GO gblDiagCol = rowMap_.getGlobalElement (lclRow);
154  const LO lclDiagCol = colMap_.getLocalElement (gblDiagCol);
155  // If it's not in the column Map, then there's no diagonal entry.
156  if (lclDiagCol != LOT::invalid ()) {
157  // TODO (mfh 25 Apr 2018) Use findRelOffset to optimize for
158  // the sorted case, but note that it requires operator[].
159  bool foundDiag = false; // don't count duplicates
160 
161  if (ignoreMapsForTriangularStructure_) {
162  for (LO k = 0; k < numEnt && ! foundDiag; ++k) {
163  const LO lclCol = G_row(k);
164  if (lclCol == lclDiagCol) {
165  foundDiag = true;
166  }
167  }
168  // mfh 30 Apr 2018: See GitHub Issue #2658. Per
169  // current Tpetra::CrsGraph::computeLocalConstants
170  // behavior, assume that local column indices are
171  // sorted in each row.
172  if (numEnt > LO (0)) {
173  const LO smallestLclCol = G_row(0);
174  const LO largestLclCol = G_row(numEnt-1); // could be same
175 
176  if (smallestLclCol < lclRow) {
177  result.couldBeUpperTriangular = false;
178  }
179  if (lclRow < largestLclCol) {
180  result.couldBeLowerTriangular = false;
181  }
182  }
183  }
184  else {
185  for (LO k = 0; k < numEnt &&
186  ((! foundDiag) ||
187  result.couldBeLowerTriangular ||
188  result.couldBeUpperTriangular);
189  ++k) {
190  const LO lclCol = G_row(k);
191  if (lclCol == lclDiagCol) {
192  foundDiag = true;
193  }
194  else {
195  const GO gblCol = colMap_.getGlobalElement (lclCol);
196  if (gblCol < gblDiagCol) {
197  result.couldBeUpperTriangular = false;
198  }
199  if (gblDiagCol < gblCol) {
200  result.couldBeLowerTriangular = false;
201  }
202  }
203  } // for each entry in lclRow
204  } // if-else ignoreMapsForTriangularStructure
205 
206  if (foundDiag) {
207  ++(result.diagCount);
208  }
209  }
210  }
211  }
212 
213  private:
214  LocalGraphType G_;
215  LocalMapType rowMap_;
216  LocalMapType colMap_;
217  bool ignoreMapsForTriangularStructure_;
218  };
219 
220 } // namespace Impl
221 
239 template<class LocalGraphType, class LocalMapType>
240 LocalTriangularStructureResult<typename LocalMapType::local_ordinal_type>
241 determineLocalTriangularStructure (const LocalGraphType& G,
242  const LocalMapType& rowMap,
243  const LocalMapType& colMap,
244  const bool ignoreMapsForTriangularStructure)
245 {
246  using LO = typename LocalMapType::local_ordinal_type;
247  using execution_space = typename LocalGraphType::device_type::execution_space;
248  using range_type = Kokkos::RangePolicy<execution_space, LO>;
249  using functor_type =
251 
252  LocalTriangularStructureResult<LO> result {0, 0, true, true};
253  Kokkos::parallel_reduce ("Tpetra::Details::determineLocalTriangularStructure",
254  range_type (0, G.numRows ()),
255  functor_type (G, rowMap, colMap,
256  ignoreMapsForTriangularStructure),
257  result);
258  return result;
259 }
260 
261 } // namespace Details
262 } // namespace Tpetra
263 
264 #endif // TPETRA_DETAILS_DETERMINELOCALTRIANGULARSTRUCTURE_HPP
Import KokkosSparse::OrdinalTraits, a traits class for &quot;invalid&quot; (flag) values of integer types...
bool couldBeLowerTriangular
Whether the graph is locally structurally lower triangular.
KOKKOS_INLINE_FUNCTION void init(result_type &dst) const
Set the initial value of the reduction result.
KOKKOS_INLINE_FUNCTION void operator()(const typename LocalMapType::local_ordinal_type lclRow, result_type &result) const
Reduction operator: result is (diagonal count, error count).
DetermineLocalTriangularStructure(const LocalGraphType &G, const LocalMapType &rowMap, const LocalMapType &colMap, const bool ignoreMapsForTriangularStructure)
Constructor.
LocalTriangularStructureResult< typename LocalMapType::local_ordinal_type > determineLocalTriangularStructure(const LocalGraphType &G, const LocalMapType &rowMap, const LocalMapType &colMap, const bool ignoreMapsForTriangularStructure)
Count the local number of diagonal entries in a local sparse graph, and determine whether the local p...
bool couldBeUpperTriangular
Whether the graph is locally structurally upper triangular.
Implementation of Tpetra::Details::determineLocalTriangularStructure (which see below).