Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
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
57namespace Tpetra {
58namespace Details {
59
63template<class LO>
74
75namespace 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
107 const LocalMapType& rowMap,
108 const LocalMapType& colMap,
110 G_ (G),
111 rowMap_ (rowMap),
112 colMap_ (colMap),
113 ignoreMapsForTriangularStructure_ (ignoreMapsForTriangularStructure)
114 {}
115
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
126 join (result_type& dst,
127 const 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
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:
215 LocalMapType rowMap_;
216 LocalMapType colMap_;
217 bool ignoreMapsForTriangularStructure_;
218 };
219
220} // namespace Impl
221
239template<class LocalGraphType, class LocalMapType>
242 const LocalMapType& rowMap,
243 const LocalMapType& colMap,
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
253 Kokkos::parallel_reduce ("Tpetra::Details::determineLocalTriangularStructure",
254 range_type (0, G.numRows ()),
255 functor_type (G, rowMap, colMap,
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 "invalid" (flag) values of integer types,...
Struct that holds views of the contents of a CrsMatrix.
Implementation of Tpetra::Details::determineLocalTriangularStructure (which see below).
KOKKOS_INLINE_FUNCTION void operator()(const typename LocalMapType::local_ordinal_type lclRow, result_type &result) const
Reduction operator: result is (diagonal count, error count).
KOKKOS_INLINE_FUNCTION void init(result_type &dst) const
Set the initial value of the reduction result.
DetermineLocalTriangularStructure(const LocalGraphType &G, const LocalMapType &rowMap, const LocalMapType &colMap, const bool ignoreMapsForTriangularStructure)
Constructor.
Implementation details of Tpetra.
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...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
bool couldBeUpperTriangular
Whether the graph is locally structurally upper triangular.
bool couldBeLowerTriangular
Whether the graph is locally structurally lower triangular.