Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_LocalCrsMatrixOperator_def.hpp
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// ************************************************************************
38// @HEADER
39
40#ifndef TPETRA_LOCALCRSMATRIXOPERATOR_DEF_HPP
41#define TPETRA_LOCALCRSMATRIXOPERATOR_DEF_HPP
42
43#include "Tpetra_LocalOperator.hpp"
45#include "KokkosSparse.hpp"
46#include "Teuchos_TestForException.hpp"
47#include "Teuchos_OrdinalTraits.hpp"
48
49namespace Tpetra {
50
51template<class MultiVectorScalar, class MatrixScalar, class Device>
52LocalCrsMatrixOperator<MultiVectorScalar, MatrixScalar, Device>::
53LocalCrsMatrixOperator (const std::shared_ptr<local_matrix_device_type>& A)
54 : A_ (A), have_A_cusparse(false)
55{
56 const char tfecfFuncName[] = "LocalCrsMatrixOperator: ";
57 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
58 (A_.get () == nullptr, std::invalid_argument,
59 "Input matrix A is null.");
60}
61
62template<class MultiVectorScalar, class MatrixScalar, class Device>
63LocalCrsMatrixOperator<MultiVectorScalar, MatrixScalar, Device>::
64LocalCrsMatrixOperator (const std::shared_ptr<local_matrix_device_type>& A, const ordinal_view_type& A_ordinal_rowptrs) :
65 A_ (A),
66 A_cusparse("LocalCrsMatrixOperator_cuSPARSE", A->numRows(), A->numCols(), A->nnz(),
67 A->values, A_ordinal_rowptrs, A->graph.entries),
68 have_A_cusparse(true)
69{
70 const char tfecfFuncName[] = "LocalCrsMatrixOperator: ";
71 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
72 (A_.get () == nullptr, std::invalid_argument,
73 "Input matrix A is null.");
74}
75
76template<class MultiVectorScalar, class MatrixScalar, class Device>
77bool
78LocalCrsMatrixOperator<MultiVectorScalar, MatrixScalar, Device>::
79hasTransposeApply () const
80{
81 return true;
82}
83
84template<class MultiVectorScalar, class MatrixScalar, class Device>
85void
86LocalCrsMatrixOperator<MultiVectorScalar, MatrixScalar, Device>::
87apply (Kokkos::View<const mv_scalar_type**, array_layout,
88 device_type, Kokkos::MemoryTraits<Kokkos::Unmanaged> > X,
89 Kokkos::View<mv_scalar_type**, array_layout,
90 device_type, Kokkos::MemoryTraits<Kokkos::Unmanaged> > Y,
91 const Teuchos::ETransp mode,
92 const mv_scalar_type alpha,
93 const mv_scalar_type beta) const
94{
95 const bool conjugate = (mode == Teuchos::CONJ_TRANS);
96 const bool transpose = (mode != Teuchos::NO_TRANS);
97
98#ifdef HAVE_TPETRA_DEBUG
99 const char tfecfFuncName[] = "apply: ";
100
101 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
102 (X.extent (1) != Y.extent (1), std::runtime_error,
103 "X.extent(1) = " << X.extent (1) << " != Y.extent(1) = "
104 << Y.extent (1) << ".");
105 // If the two pointers are NULL, then they don't alias one
106 // another, even though they are equal.
107 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
108 (X.data () == Y.data () && X.data () != nullptr,
109 std::runtime_error, "X and Y may not alias one another.");
110#endif // HAVE_TPETRA_DEBUG
111
112 const auto op = transpose ?
113 (conjugate ? KokkosSparse::ConjugateTranspose :
114 KokkosSparse::Transpose) : KokkosSparse::NoTranspose;
115 //Currently KK has no cusparse wrapper for rank-2 (SpMM)
116 //TODO: whent that is supported, use A_cusparse for that case also
117 if(X.extent(1) == size_t(1) && have_A_cusparse)
118 {
119 KokkosSparse::spmv (op, alpha, A_cusparse, Kokkos::subview(X, Kokkos::ALL(), 0),
120 beta, Kokkos::subview(Y, Kokkos::ALL(), 0));
121 }
122 else
123 {
124 KokkosSparse::spmv (op, alpha, *A_, X, beta, Y);
125 }
126}
127
130template<class MultiVectorScalar, class MatrixScalar, class Device>
131void
134 Kokkos::View<const mv_scalar_type**, array_layout,
135 device_type, Kokkos::MemoryTraits<Kokkos::Unmanaged> > X,
136 Kokkos::View<mv_scalar_type**, array_layout,
137 device_type, Kokkos::MemoryTraits<Kokkos::Unmanaged> > Y,
138 const Teuchos::ETransp mode,
139 const mv_scalar_type alpha,
140 const mv_scalar_type beta) const
141{
142 const bool conjugate = (mode == Teuchos::CONJ_TRANS);
143 const bool transpose = (mode != Teuchos::NO_TRANS);
144
145#ifdef HAVE_TPETRA_DEBUG
146 const char tfecfFuncName[] = "applyLoadBalanced: ";
147
149 (X.extent (1) != Y.extent (1), std::runtime_error,
150 "X.extent(1) = " << X.extent (1) << " != Y.extent(1) = "
151 << Y.extent (1) << ".");
152 // If the two pointers are NULL, then they don't alias one
153 // another, even though they are equal.
155 (X.data () == Y.data () && X.data () != nullptr,
156 std::runtime_error, "X and Y may not alias one another.");
157#endif // HAVE_TPETRA_DEBUG
158
159 const auto op = transpose ?
160 (conjugate ? KokkosSparse::ConjugateTranspose :
161 KokkosSparse::Transpose) : KokkosSparse::NoTranspose;
162 //Select the merge path algorithm (used if available, otherwise has no effect)
163 //TODO BMK: If/when KokkosKernels gets its own SPMV implementation for imbalanced rows,
164 //call that here or select it using Controls.
165 //Ideally it supports multivectors from the beginning.
166 if((Details::Behavior::useMergePathMultiVector() || X.extent(1) == size_t(1)) && have_A_cusparse)
167 {
168 KokkosKernels::Experimental::Controls controls;
169 controls.setParameter("algorithm", "merge");
170 //Apply on one column at a time (must be rank-1)
171 for(size_t vec = 0; vec < X.extent(1); vec++)
172 {
173 KokkosSparse::spmv (controls, op,
174 alpha, A_cusparse, Kokkos::subview(X, Kokkos::ALL(), vec),
175 beta, Kokkos::subview(Y, Kokkos::ALL(), vec));
176 }
177 }
178 else
179 {
180 //Just run multivector version of spmv (no controls, and no cusparse support)
181 KokkosSparse::spmv (op, alpha, *A_, X, beta, Y);
182 }
183}
184
185template<class MultiVectorScalar, class MatrixScalar, class Device>
186const typename LocalCrsMatrixOperator<MultiVectorScalar, MatrixScalar, Device>::local_matrix_device_type&
189{
190 return *A_;
191}
192
193} // namespace Tpetra
194
195//
196// Explicit instantiation macro
197//
198// Must be expanded from within the Tpetra namespace!
199//
200
201// We only explicitly instantiate for MultiVectorScalar ==
202// MatrixScalar, which is what CrsMatrix needs.
203
204#define TPETRA_LOCALCRSMATRIXOPERATOR_INSTANT(SC,NT) \
205 template class LocalCrsMatrixOperator< SC, SC, NT::device_type >;
206
207// If we want mixed versions, we use this macro.
208
209#define TPETRA_LOCALCRSMATRIXOPERATOR_MIXED_INSTANT(SC,MATSC,LO,GO,NT) \
210 template class LocalCrsMatrixOperator< SC, MATSC, NT::device_type >;
211
212#endif // TPETRA_LOCALCRSMATRIXOPERATOR_DEF_HPP
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra's behavior.
Struct that holds views of the contents of a CrsMatrix.
Abstract interface for local operators (e.g., matrices and preconditioners).
Namespace Tpetra contains the class and methods constituting the Tpetra library.