Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tpetra_Apply_Helpers.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_APPLY_HELPERS_HPP
41 #define TPETRA_APPLY_HELPERS_HPP
42 #include <type_traits>
43 #include "Teuchos_ScalarTraits.hpp"
44 #include "Tpetra_CrsMatrix.hpp"
46 
47 namespace Tpetra {
49 
50 
73  template <class MatrixArray, class MultiVectorArray>
74  void batchedApply(const MatrixArray &Matrices,
75  const typename std::remove_pointer<typename MultiVectorArray::value_type>::type &X,
76  MultiVectorArray &Y,
77  typename std::remove_pointer<typename MatrixArray::value_type>::type::scalar_type alpha = Teuchos::ScalarTraits< typename std::remove_pointer<typename MatrixArray::value_type>::type::scalar_type>::one(),
78  typename std::remove_pointer<typename MatrixArray::value_type>::type::scalar_type beta = Teuchos::ScalarTraits< typename std::remove_pointer<typename MatrixArray::value_type>::type::scalar_type>::zero(),
79  Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::null) {
80  Tpetra::Details::ProfilingRegion regionTotal ("Tpetra::batchedApply");
81 
82  using size_type = typename MatrixArray::size_type;
83  using matrix_type = typename std::remove_pointer<typename MatrixArray::value_type>::type;
84  using map_type = typename matrix_type::map_type;
85  using import_type = typename matrix_type::import_type;
86  using export_type = typename matrix_type::export_type;
87  using MV = typename matrix_type::MV;
88  using scalar_type = typename matrix_type::scalar_type;
89 
90  using Teuchos::RCP;
91  using Teuchos::rcp_const_cast;
92 
93  const scalar_type ONE = Teuchos::ScalarTraits<scalar_type>::one();
94  const scalar_type ZERO = Teuchos::ScalarTraits<scalar_type>::zero();
95 
96  size_type N = Matrices.size();
97  if(N==0) return;
98  int numRanks = X.getMap()->getComm()->getSize();
99 
100  // If X is aliased to any Y but the last one, throw
101  for(size_type i=0; i<N-1; i++) {
102  TEUCHOS_TEST_FOR_EXCEPTION( &X == Y[i], std::runtime_error, "Tpetra::batchedApply(): X cannot be aliased to any Y except the final one.");
103  }
104 
105  /* Checks for whether the reduced communication path can be used */
106  RCP<const map_type> compare_colMap = Matrices[0]->getColMap();
107  RCP<const import_type> importer = Matrices[0]->getGraph()->getImporter();
108 
109 
110  bool can_batch, check_maps;
111  if(params.is_null() || !params->isParameter("can batch")) {
112  can_batch = (importer.is_null() || N==1) ? false :true;
113  check_maps = true;
114  }
115  else {
116  can_batch = (!params->get<bool>("can batch") || importer.is_null() || N==1) ? false : true;
117  check_maps = false;
118  }
119  // We can't batch with replicated Y's
120  if(numRanks > 1) {
121  for(size_type i=0; i<N && can_batch; i++) {
122  if(!Y[i]->isDistributed())
123  can_batch = false;
124  }
125  }
126 
127  // Do the domain/column maps all match?
128  for(size_type i=1; i<N && check_maps && can_batch; i++) {
129  if(!Matrices[i]->getColMap()->isSameAs(*compare_colMap)) {
130  can_batch=false;
131  }
132  }
133 
134  if(can_batch) {
135  /* Batching path: Guarantees an existing importer and N>1 */
136 
137  // Special case for alpha = 0
138  if (alpha == ZERO) {
139  if (beta == ZERO) {
140  for(size_type i=0; i<N; i++) Y[i]->putScalar(ZERO);
141  } else if (beta != ONE) {
142  for(size_type i=0; i<N; i++) Y[i]->scale(beta);
143  }
144  if(!params.is_null()) params->set("can batch",true);
145  return;
146  }
147 
148  const bool Y_is_overwritten = (beta == ZERO);
149 
150  // Start by importing X to Matrices[0]'s temporary
151  RCP<const MV> X_colMap;
152  {
153  Tpetra::Details::ProfilingRegion regionImport ("Tpetra::batchedApply: Import");
154  RCP<MV> X_colMapNonConst = Matrices[0]->getColumnMapMultiVector(X);
155 
156  // Import from the domain Map MV to the column Map MV.
157  X_colMapNonConst->doImport(X, *importer, INSERT);
158  X_colMap = rcp_const_cast<const MV>(X_colMapNonConst);
159  }
160 
161  for(size_type i=0; i<N; i++) {
162  RCP<const export_type> exporter = Matrices[i]->getGraph()->getExporter();
163 
164  // Temporary MV for doExport (if needed),
165  RCP<MV> Y_rowMap = Matrices[i]->getRowMapMultiVector(*Y[i]);
166  if (!exporter.is_null()) {
167  Matrices[i]->localApply(*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha, beta);
168  {
169  Tpetra::Details::ProfilingRegion regionExport ("Tpetra::batchedApply: Export");
170  if (Y_is_overwritten) {
171  Y[i]->putScalar(ZERO);
172  }
173  else {
174  Y[i]->scale (beta);
175  }
176  Y[i]->doExport(*Y_rowMap, *exporter, ADD);
177  }
178  }
179  else { // Don't do an Export: row Map and range Map are the same.
180  // Check for aliasing
181  if (! Y[i]->isConstantStride() || X_colMap.getRawPtr() == Y[i]) {
182  Y_rowMap = Matrices[i]->getRowMapMultiVector(*Y[i], true);
183  if (beta != ZERO) {
184  Tpetra::deep_copy (*Y_rowMap, *Y[i]);
185  }
186 
187  Matrices[i]->localApply(*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha, beta);
188  Tpetra::deep_copy(*Y[i], *Y_rowMap);
189  }
190  else {
191  Matrices[i]->localApply(*X_colMap, *Y[i], Teuchos::NO_TRANS, alpha, beta);
192  }
193  }
194  }
195  if(!params.is_null()) params->set("can batch",true);
196  }
197  else {
198  /* Non-batching path */
199  for(size_type i=0; i<N; i++) {
200  Matrices[i]->apply(X,*Y[i],Teuchos::NO_TRANS, alpha, beta);
201  }
202  if(!params.is_null()) params->set("can batch",false);
203  }
204  }
206 
207 }// namespace Tpetra
208 
209 #endif // TPETRA_APPLY_HELPERS_HPP
210 
Declaration of Tpetra::Details::Profiling, a scope guard for Kokkos Profiling.
void batchedApply(const MatrixArray &Matrices, const typename std::remove_pointer< typename MultiVectorArray::value_type >::type &X, MultiVectorArray &Y, typename std::remove_pointer< typename MatrixArray::value_type >::type::scalar_type alpha=Teuchos::ScalarTraits< typename std::remove_pointer< typename MatrixArray::value_type >::type::scalar_type >::one(), typename std::remove_pointer< typename MatrixArray::value_type >::type::scalar_type beta=Teuchos::ScalarTraits< typename std::remove_pointer< typename MatrixArray::value_type >::type::scalar_type >::zero(), Teuchos::RCP< Teuchos::ParameterList > params=Teuchos::null)
Does multiply matrix apply() calls with a single X vector.
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.
Insert new values that don&#39;t currently exist.
Sum new values into existing values.
Replace old values with zero.