Tpetra parallel linear algebra  Version of the Day
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"
61 #include "Tpetra_Vector_decl.hpp"
62 #include <type_traits>
63 
64 namespace Tpetra {
65 namespace Details {
66 
79 template<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  {
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 
144 private:
146  DiagType D_;
148  LocalMapType rowMap_;
150  LocalMapType colMap_;
152  CrsMatrixType A_;
153 };
154 
155 
178 template<class DiagType,
179  class LocalMapType,
180  class CrsMatrixType>
181 static typename LocalMapType::local_ordinal_type
182 getDiagCopyWithoutOffsets (const DiagType& D,
183  const LocalMapType& rowMap,
184  const LocalMapType& colMap,
185  const CrsMatrixType& A)
186 {
187  static_assert (Kokkos::Impl::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<DiagType, typename DiagType::non_const_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.dimension_0 ());
206  LO errCount = 0;
207  Kokkos::parallel_reduce (policy_type (0, numRows), functor, errCount);
208  return errCount;
209 }
210 
247 template<class SC, class LO, class GO, class NT>
248 LO
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
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Traits class for "invalid" (flag) values of integer types that Tpetra uses as local ordinals or globa...
Functor that implements much of the one-argument overload of Tpetra::CrsMatrix::getLocalDiagCopy, for the case where the matrix is fill complete.
LO value_type
The result of the reduction; number of errors.
Declaration of the Tpetra::Vector class.
KOKKOS_FUNCTION void operator()(const LO &lclRowInd, value_type &errCount) const
Operator for Kokkos::parallel_for.
Implementation details of Tpetra.
Traits class for "invalid" (flag) values of integer types that Tpetra uses as local ordinals or globa...
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...
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&#39;s diagonal entries into a Tpetra::V...
Replace old values with zero.
A distributed dense vector.
CrsMatrixGetDiagCopyFunctor(const DiagType &D, const LocalMapType &rowMap, const LocalMapType &colMap, const CrsMatrixType &A)
Constructor.