Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
58 namespace Tpetra {
59 namespace Details {
60 
61 template<class RV, class XMV>
62 void
63 lclDot (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::Impl::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::Impl::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  Kokkos::deep_copy (theDots, zero);
149  }
150  else { // lclNumRows != 0
151  if (constantStrideX && constantStrideY) {
152  if (X.extent (1) == 1) {
153  typename RV::non_const_value_type result =
154  KokkosBlas::dot (subview (X, ALL (), 0), subview (Y, ALL (), 0));
155  Kokkos::deep_copy (theDots, result);
156  }
157  else {
158  KokkosBlas::dot (theDots, X, Y);
159  }
160  }
161  else { // not constant stride
162  // NOTE (mfh 15 Jul 2014) This does a kernel launch for
163  // every column. It might be better to have a kernel that
164  // does the work all at once. On the other hand, we don't
165  // prioritize performance of MultiVector views of
166  // noncontiguous columns.
167  for (size_t k = 0; k < numVecs; ++k) {
168  const size_t X_col = constantStrideX ? k : whichVecsX[k];
169  const size_t Y_col = constantStrideY ? k : whichVecsY[k];
170  KokkosBlas::dot (subview (theDots, k), subview (X, ALL (), X_col),
171  subview (Y, ALL (), Y_col));
172  } // for each column
173  } // constantStride
174  } // lclNumRows != 0
175 }
176 
177 } // namespace Details
178 } // namespace Tpetra
179 
180 #endif // TPETRA_DETAILS_LCLDOT_HPP
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
Copy the contents of the MultiVector src into dst.