Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_Details_getDiagCopyWithoutOffsets_decl.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_GETDIAGCOPYWITHOUTOFFSETS_DECL_HPP
45#define TPETRA_DETAILS_GETDIAGCOPYWITHOUTOFFSETS_DECL_HPP
46
55
56#include "TpetraCore_config.h"
57#include "Kokkos_Core.hpp"
58#include "Kokkos_ArithTraits.hpp"
60#include "Tpetra_RowMatrix_decl.hpp"
62#include <type_traits>
63
64namespace Tpetra {
65namespace Details {
66
79template<class DiagType,
80 class LocalMapType,
81 class CrsMatrixType>
83 typedef typename LocalMapType::local_ordinal_type LO; // local ordinal type
84 typedef typename LocalMapType::global_ordinal_type GO; // global ordinal type
85 typedef typename CrsMatrixType::device_type device_type;
86 typedef typename CrsMatrixType::value_type scalar_type;
87 typedef typename CrsMatrixType::size_type offset_type;
88
90 typedef LO value_type;
91
99 CrsMatrixGetDiagCopyFunctor (const DiagType& D,
100 const LocalMapType& rowMap,
101 const LocalMapType& colMap,
102 const CrsMatrixType& A) :
103 D_ (D), rowMap_ (rowMap), colMap_ (colMap), A_ (A)
104 {}
105
109 KOKKOS_FUNCTION void
110 operator () (const LO& lclRowInd, value_type& errCount) const
111 {
112 const LO INV = Tpetra::Details::OrdinalTraits<LO>::invalid ();
113 const scalar_type ZERO =
114 Kokkos::Details::ArithTraits<scalar_type>::zero ();
115
116 // If the row lacks a stored diagonal entry, then its value is zero.
117 D_(lclRowInd) = ZERO;
118 const GO gblInd = rowMap_.getGlobalElement (lclRowInd);
119 const LO lclColInd = colMap_.getLocalElement (gblInd);
120
121 if (lclColInd != INV) {
122 auto curRow = A_.rowConst (lclRowInd);
123
124 // FIXME (mfh 12 May 2016) Use binary search when the row is
125 // long enough. findRelOffset currently lives in KokkosKernels
126 // (in tpetra/kernels/src/Kokkos_Sparse_findRelOffset.hpp).
127 LO offset = 0;
128 const LO numEnt = curRow.length;
129 for ( ; offset < numEnt; ++offset) {
130 if (curRow.colidx(offset) == lclColInd) {
131 break;
132 }
133 }
134
135 if (offset == numEnt) {
136 ++errCount;
137 }
138 else {
139 D_(lclRowInd) = curRow.value(offset);
140 }
141 }
142 }
143
144private:
146 DiagType D_;
148 LocalMapType rowMap_;
150 LocalMapType colMap_;
152 CrsMatrixType A_;
153};
154
155
178template<class DiagType,
179 class LocalMapType,
180 class CrsMatrixType>
181static typename LocalMapType::local_ordinal_type
182getDiagCopyWithoutOffsets (const DiagType& D,
183 const LocalMapType& rowMap,
184 const LocalMapType& colMap,
185 const CrsMatrixType& A)
186{
187 static_assert (Kokkos::is_view<DiagType>::value,
188 "DiagType must be a Kokkos::View.");
189 static_assert (static_cast<int> (DiagType::rank) == 1,
190 "DiagType must be a 1-D Kokkos::View.");
191 static_assert (std::is_same<typename DiagType::value_type, typename DiagType::non_const_value_type>::value,
192 "DiagType must be a nonconst Kokkos::View.");
193 typedef typename LocalMapType::local_ordinal_type LO;
194 typedef typename CrsMatrixType::device_type device_type;
195 typedef typename device_type::execution_space execution_space;
196 typedef Kokkos::RangePolicy<execution_space, LO> policy_type;
197
198 typedef Kokkos::View<typename DiagType::non_const_value_type*,
199 typename DiagType::array_layout,
200 typename DiagType::device_type,
201 Kokkos::MemoryUnmanaged> diag_type;
202 diag_type D_out = D;
204 functor (D_out, rowMap, colMap, A);
205 const LO numRows = static_cast<LO> (D.extent (0));
206 LO errCount = 0;
207 Kokkos::parallel_reduce (policy_type (0, numRows), functor, errCount);
208 return errCount;
209}
210
247template<class SC, class LO, class GO, class NT>
248LO
250 const ::Tpetra::RowMatrix<SC, LO, GO, NT>& A,
251 const bool debug =
252#ifdef HAVE_TPETRA_DEBUG
253 true);
254#else // ! HAVE_TPETRA_DEBUG
255 false);
256#endif // HAVE_TPETRA_DEBUG
257
258} // namespace Details
259} // namespace Tpetra
260
261#endif // TPETRA_DETAILS_GETDIAGCOPYWITHOUTOFFSETS_DECL_HPP
Import KokkosSparse::OrdinalTraits, a traits class for "invalid" (flag) values of integer types,...
Declaration of the Tpetra::Vector class.
A distributed dense vector.
Implementation details of Tpetra.
LO getLocalDiagCopyWithoutOffsetsNotFillComplete(::Tpetra::Vector< SC, LO, GO, NT > &diag, const ::Tpetra::RowMatrix< SC, LO, GO, NT > &A, const bool debug=false)
Given a locally indexed, global sparse matrix, extract the matrix's diagonal entries into a Tpetra::V...
static LocalMapType::local_ordinal_type getDiagCopyWithoutOffsets(const DiagType &D, const LocalMapType &rowMap, const LocalMapType &colMap, const CrsMatrixType &A)
Given a locally indexed, local sparse matrix, and corresponding local row and column Maps,...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
@ ZERO
Replace old values with zero.
Functor that implements much of the one-argument overload of Tpetra::CrsMatrix::getLocalDiagCopy,...
CrsMatrixGetDiagCopyFunctor(const DiagType &D, const LocalMapType &rowMap, const LocalMapType &colMap, const CrsMatrixType &A)
Constructor.
KOKKOS_FUNCTION void operator()(const LO &lclRowInd, value_type &errCount) const
Operator for Kokkos::parallel_for.