Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_Details_lclDot.hpp
Go to the documentation of this file.
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// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41
42#ifndef TPETRA_DETAILS_LCLDOT_HPP
43#define TPETRA_DETAILS_LCLDOT_HPP
44
51
52#include "Kokkos_DualView.hpp"
53#include "Kokkos_ArithTraits.hpp"
54#include "KokkosBlas1_dot.hpp"
55#include "Teuchos_ArrayView.hpp"
56#include "Teuchos_TestForException.hpp"
57
58namespace Tpetra {
59namespace Details {
60
61template<class RV, class XMV>
62void
63lclDot (const RV& dotsOut,
64 const XMV& X_lcl,
65 const XMV& Y_lcl,
66 const size_t lclNumRows,
67 const size_t numVecs,
68 const size_t whichVecsX[],
69 const size_t whichVecsY[],
70 const bool constantStrideX,
71 const bool constantStrideY)
72{
73 using Kokkos::ALL;
74 using Kokkos::subview;
75 typedef typename RV::non_const_value_type dot_type;
76#ifdef HAVE_TPETRA_DEBUG
77 const char prefix[] = "Tpetra::MultiVector::lclDotImpl: ";
78#endif // HAVE_TPETRA_DEBUG
79
80 static_assert (Kokkos::is_view<RV>::value,
81 "Tpetra::MultiVector::lclDotImpl: "
82 "The first argument dotsOut is not a Kokkos::View.");
83 static_assert (RV::rank == 1, "Tpetra::MultiVector::lclDotImpl: "
84 "The first argument dotsOut must have rank 1.");
85 static_assert (Kokkos::is_view<XMV>::value,
86 "Tpetra::MultiVector::lclDotImpl: The type of the 2nd and "
87 "3rd arguments (X_lcl and Y_lcl) is not a Kokkos::View.");
88 static_assert (XMV::rank == 2, "Tpetra::MultiVector::lclDotImpl: "
89 "X_lcl and Y_lcl must have rank 2.");
90
91 // In case the input dimensions don't match, make sure that we
92 // don't overwrite memory that doesn't belong to us, by using
93 // subset views with the minimum dimensions over all input.
94 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
95 const std::pair<size_t, size_t> colRng (0, numVecs);
96 RV theDots = subview (dotsOut, colRng);
97 XMV X = subview (X_lcl, rowRng, ALL ());
98 XMV Y = subview (Y_lcl, rowRng, ALL ());
99
100#ifdef HAVE_TPETRA_DEBUG
101 if (lclNumRows != 0) {
102 TEUCHOS_TEST_FOR_EXCEPTION
103 (X.extent (0) != lclNumRows, std::logic_error, prefix <<
104 "X.extent(0) = " << X.extent (0) << " != lclNumRows "
105 "= " << lclNumRows << ". "
106 "Please report this bug to the Tpetra developers.");
107 TEUCHOS_TEST_FOR_EXCEPTION
108 (Y.extent (0) != lclNumRows, std::logic_error, prefix <<
109 "Y.extent(0) = " << Y.extent (0) << " != lclNumRows "
110 "= " << lclNumRows << ". "
111 "Please report this bug to the Tpetra developers.");
112 // If a MultiVector is constant stride, then numVecs should
113 // equal its View's number of columns. Otherwise, numVecs
114 // should be less than its View's number of columns.
115 TEUCHOS_TEST_FOR_EXCEPTION
116 (constantStrideX &&
117 (X.extent (0) != lclNumRows || X.extent (1) != numVecs),
118 std::logic_error, prefix << "X is " << X.extent (0) << " x " <<
119 X.extent (1) << " (constant stride), which differs from the "
120 "local dimensions " << lclNumRows << " x " << numVecs << ". "
121 "Please report this bug to the Tpetra developers.");
122 TEUCHOS_TEST_FOR_EXCEPTION
123 (! constantStrideX &&
124 (X.extent (0) != lclNumRows || X.extent (1) < numVecs),
125 std::logic_error, prefix << "X is " << X.extent (0) << " x " <<
126 X.extent (1) << " (NOT constant stride), but the local "
127 "dimensions are " << lclNumRows << " x " << numVecs << ". "
128 "Please report this bug to the Tpetra developers.");
129 TEUCHOS_TEST_FOR_EXCEPTION
130 (constantStrideY &&
131 (Y.extent (0) != lclNumRows || Y.extent (1) != numVecs),
132 std::logic_error, prefix << "Y is " << Y.extent (0) << " x " <<
133 Y.extent (1) << " (constant stride), which differs from the "
134 "local dimensions " << lclNumRows << " x " << numVecs << ". "
135 "Please report this bug to the Tpetra developers.");
136 TEUCHOS_TEST_FOR_EXCEPTION
137 (! constantStrideY &&
138 (Y.extent (0) != lclNumRows || Y.extent (1) < numVecs),
139 std::logic_error, prefix << "Y is " << Y.extent (0) << " x " <<
140 Y.extent (1) << " (NOT constant stride), but the local "
141 "dimensions are " << lclNumRows << " x " << numVecs << ". "
142 "Please report this bug to the Tpetra developers.");
143 }
144#endif // HAVE_TPETRA_DEBUG
145
146 if (lclNumRows == 0) {
147 const dot_type zero = Kokkos::Details::ArithTraits<dot_type>::zero ();
148 // DEEP_COPY REVIEW - NOT TESTED
149 Kokkos::deep_copy (theDots, zero);
150 }
151 else { // lclNumRows != 0
152 if (constantStrideX && constantStrideY) {
153 if (X.extent (1) == 1) {
154 typename RV::non_const_value_type result =
155 KokkosBlas::dot (subview (X, ALL (), 0), subview (Y, ALL (), 0));
156 // DEEP_COPY REVIEW - NOT TESTED
157 Kokkos::deep_copy (theDots, result);
158 }
159 else {
160 KokkosBlas::dot (theDots, X, Y);
161 }
162 }
163 else { // not constant stride
164 // NOTE (mfh 15 Jul 2014) This does a kernel launch for
165 // every column. It might be better to have a kernel that
166 // does the work all at once. On the other hand, we don't
167 // prioritize performance of MultiVector views of
168 // noncontiguous columns.
169 for (size_t k = 0; k < numVecs; ++k) {
170 const size_t X_col = constantStrideX ? k : whichVecsX[k];
171 const size_t Y_col = constantStrideY ? k : whichVecsY[k];
172 KokkosBlas::dot (subview (theDots, k), subview (X, ALL (), X_col),
173 subview (Y, ALL (), Y_col));
174 } // for each column
175 } // constantStride
176 } // lclNumRows != 0
177}
178
179} // namespace Details
180} // namespace Tpetra
181
182#endif // TPETRA_DETAILS_LCLDOT_HPP
Implementation details of Tpetra.
Namespace Tpetra contains the class and methods constituting the Tpetra library.