Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
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
47namespace Tpetra {
49
50
72
73 template <class MatrixArray, class MultiVectorArray>
75 const typename std::remove_pointer<typename MultiVectorArray::value_type>::type &X,
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 */
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
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);
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, ZERO);
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_ASSIGN);
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) {
185 }
186
187 Matrices[i]->localApply(*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha, beta);
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.
Struct that holds views of the contents of a CrsMatrix.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
A parallel distribution of indices over processes.
One or more distributed dense vectors.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
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.
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.
@ ADD_ASSIGN
Accumulate new values into existing values (may not be supported in all classes)
@ INSERT
Insert new values that don't currently exist.
@ ZERO
Replace old values with zero.