Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Details_MultiVectorLocalGatherScatter.hpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
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
43#ifndef IFPACK2_DETAILS_MULTIVECTORLOCALGATHERSCATTER_HPP
44#define IFPACK2_DETAILS_MULTIVECTORLOCALGATHERSCATTER_HPP
45
49
50#include "Tpetra_MultiVector.hpp"
51#include "Tpetra_Map.hpp"
52
53namespace Ifpack2 {
54namespace Details {
55
84template<class MV_in, class MV_out>
86public:
87 typedef typename MV_in::scalar_type InScalar;
88 typedef typename MV_out::scalar_type OutScalar;
89 typedef typename MV_in::local_ordinal_type LO;
90 typedef typename MV_in::global_ordinal_type GO;
91 typedef typename MV_in::node_type NO;
92
93 /**************/
94 /* MV <==> MV */
95 /**************/
96 void
97 gather (MV_out& X_out,
98 const MV_in& X_in,
99 const Teuchos::ArrayView<const LO> perm) const
100 {
101 using Teuchos::ArrayRCP;
102 const size_t numRows = X_out.getLocalLength ();
103 const size_t numVecs = X_in.getNumVectors ();
104 Kokkos::fence();
105 for (size_t j = 0; j < numVecs; ++j) {
106 ArrayRCP<const InScalar> X_in_j = X_in.getData(j);
107 ArrayRCP<OutScalar> X_out_j = X_out.getDataNonConst(j);
108 for (size_t i = 0; i < numRows; ++i) {
109 const size_t i_perm = perm[i];
110 X_out_j[i] = X_in_j[i_perm];
111 }
112 }
113 X_out.modify_host();
114 }
115
116 //Gather blocks (contiguous groups of blockSize rows)
117 //X_out and X_in are point indexed, but perm uses block indices.
118 //So X_out.getLocalLength() / blockSize gives the number of blocks.
119 void
120 gatherBlock (
121 MV_out& X_out,
122 const MV_in& X_in,
123 const Teuchos::ArrayView<const LO> perm,
124 LO blockSize) const
125 {
126 using Teuchos::ArrayRCP;
127 const size_t numBlocks = X_out.getLocalLength() / blockSize;
128 const size_t numVecs = X_in.getNumVectors ();
129 Kokkos::fence();
130 for (size_t j = 0; j < numVecs; ++j) {
131 ArrayRCP<const InScalar> X_in_j = X_in.getData(j);
132 ArrayRCP<OutScalar> X_out_j = X_out.getDataNonConst(j);
133 for (size_t i = 0; i < numBlocks; ++i) {
134 const size_t i_perm = perm[i];
135 for (LO k = 0; k < blockSize; k++) {
136 X_out_j[i * blockSize + k] = X_in_j[i_perm * blockSize + k];
137 }
138 }
139 }
140 X_out.modify_host();
141 }
142
143 void
144 scatter (MV_in& X_in,
145 const MV_out& X_out,
146 const Teuchos::ArrayView<const LO> perm) const
147 {
148 using Teuchos::ArrayRCP;
149 const size_t numRows = X_out.getLocalLength();
150 const size_t numVecs = X_in.getNumVectors();
151 Kokkos::fence();
152 for (size_t j = 0; j < numVecs; ++j) {
153 ArrayRCP<InScalar> X_in_j = X_in.getDataNonConst(j);
154 ArrayRCP<const OutScalar> X_out_j = X_out.getData(j);
155 for (size_t i = 0; i < numRows; ++i) {
156 const size_t i_perm = perm[i];
157 X_in_j[i_perm] = X_out_j[i];
158 }
159 }
160 X_out.modify_host();
161 }
162
163 void
164 scatterBlock (
165 MV_in& X_in,
166 const MV_out& X_out,
167 const Teuchos::ArrayView<const LO> perm,
168 LO blockSize) const
169 {
170 using Teuchos::ArrayRCP;
171 const size_t numBlocks = X_out.getLocalLength() / blockSize;
172 const size_t numVecs = X_in.getNumVectors ();
173 Kokkos::fence();
174 for (size_t j = 0; j < numVecs; ++j) {
175 ArrayRCP<const InScalar> X_in_j = X_in.getData(j);
176 ArrayRCP<OutScalar> X_out_j = X_out.getDataNonConst(j);
177 for (size_t i = 0; i < numBlocks; ++i) {
178 const size_t i_perm = perm[i];
179 for (LO k = 0; k < blockSize; k++) {
180 X_in_j[i_perm * blockSize + k] = X_out_j[i * blockSize + k];
181 }
182 }
183 }
184 X_out.modify_host();
185 }
186
187 /******************/
188 /* View <==> View */
189 /******************/
190 template<typename InView, typename OutView>
191 void gatherViewToView(OutView X_out,
192 const InView X_in,
193 const Teuchos::ArrayView<const LO> perm) const
194 {
195 //note: j is col, i is row
196 Kokkos::fence(); // demonstrated via unit test failure
197 for(size_t j = 0; j < X_out.extent(1); ++j) {
198 for(size_t i = 0; i < X_out.extent(0); ++i) {
199 const LO i_perm = perm[i];
200 X_out(i, j) = X_in(i_perm, j);
201 }
202 }
203 }
204
205 template<typename InView, typename OutView>
206 void scatterViewToView(InView X_in,
207 const OutView X_out,
208 const Teuchos::ArrayView<const LO> perm) const
209 {
210 Kokkos::fence();
211 for(size_t j = 0; j < X_out.extent(1); ++j) {
212 for(size_t i = 0; i < X_out.extent(0); ++i) {
213 const LO i_perm = perm[i];
214 X_in(i_perm, j) = X_out(i, j);
215 }
216 }
217 }
218
219 template<typename InView, typename OutView>
220 void gatherViewToViewBlock(OutView X_out,
221 const InView X_in,
222 const Teuchos::ArrayView<const LO> perm,
223 LO blockSize) const
224 {
225 //note: j is col, i is row
226 Kokkos::fence();
227 size_t numBlocks = X_out.extent(0) / blockSize;
228 for(size_t j = 0; j < X_out.extent(1); ++j) {
229 for(size_t i = 0; i < numBlocks; ++i) {
230 const LO i_perm = perm[i];
231 for(LO k = 0; k < blockSize; k++) {
232 X_out(i * blockSize + k, j) = X_in(i_perm * blockSize + k, j);
233 }
234 }
235 }
236 }
237
238 template<typename InView, typename OutView>
239 void scatterViewToViewBlock(InView X_in,
240 const OutView X_out,
241 const Teuchos::ArrayView<const LO> perm,
242 LO blockSize) const
243 {
244 //note: j is col, i is row
245 Kokkos::fence();
246 size_t numBlocks = X_out.extent(0) / blockSize;
247 for(size_t j = 0; j < X_out.extent(1); ++j) {
248 for(size_t i = 0; i < numBlocks; ++i) {
249 const LO i_perm = perm[i];
250 for(LO k = 0; k < blockSize; k++) {
251 X_in(i_perm * blockSize + k, j) = X_out(i * blockSize + k, j);
252 }
253 }
254 }
255 }
256
257 /*******************************/
258 /* MV <==> View specialization */
259 /*******************************/
260 template<typename InView>
261 void gatherMVtoView(MV_out X_out,
262 InView X_in,
263 const Teuchos::ArrayView<const LO> perm) const
264 {
265 //note: j is col, i is row
266 Kokkos::fence();
267 size_t numRows = X_out.getLocalLength();
268 for(size_t j = 0; j < X_out.getNumVectors(); ++j) {
269 Teuchos::ArrayRCP<OutScalar> X_out_j = X_out.getDataNonConst(j);
270 for(size_t i = 0; i < numRows; ++i) {
271 const LO i_perm = perm[i];
272 X_out_j[i] = X_in(i_perm, j);
273 }
274 }
275 }
276
277 template<typename InView>
278 void scatterMVtoView(InView X_in,
279 MV_out X_out,
280 const Teuchos::ArrayView<const LO> perm) const
281 {
282 size_t numRows = X_out.getLocalLength();
283 Kokkos::fence();
284 for(size_t j = 0; j < X_in.extent(1); ++j) {
285 Teuchos::ArrayRCP<const OutScalar> X_out_j = X_out.getData(j);
286 for(size_t i = 0; i < numRows; ++i) {
287 const LO i_perm = perm[i];
288 X_in(i_perm, j) = X_out_j[i];
289 }
290 }
291 }
292
293 template<typename InView>
294 void gatherMVtoViewBlock(MV_out X_out,
295 InView X_in,
296 const Teuchos::ArrayView<const LO> perm,
297 LO blockSize) const
298 {
299 //note: j is col, i is row
300 size_t numBlocks = X_out.getLocalLength() / blockSize;
301 Kokkos::fence();
302 for(size_t j = 0; j < X_out.getNumVectors(); ++j) {
303 Teuchos::ArrayRCP<OutScalar> X_out_j = X_out.getDataNonConst(j);
304 for(size_t i = 0; i < numBlocks; ++i) {
305 const LO i_perm = perm[i];
306 for(LO k = 0; k < blockSize; k++) {
307 X_out_j[i * blockSize + k] = X_in(i_perm * blockSize + k, j);
308 }
309 }
310 }
311 }
312
313 template<typename InView>
314 void scatterMVtoViewBlock(InView X_in,
315 MV_out X_out,
316 const Teuchos::ArrayView<const LO> perm,
317 LO blockSize) const
318 {
319 size_t numBlocks = X_out.getLocalLength() / blockSize;
320 Kokkos::fence();
321 for(size_t j = 0; j < X_in.extent(1); ++j) {
322 Teuchos::ArrayRCP<const OutScalar> X_out_j = X_out.getData(j);
323 for(size_t i = 0; i < numBlocks; ++i) {
324 const LO i_perm = perm[i];
325 for(LO k = 0; k < blockSize; k++) {
326 X_in(i_perm * blockSize + k, j) = X_out_j[i * blockSize + k];
327 }
328 }
329 }
330 }
331};
332
333} // namespace Details
334} // namespace Ifpack2
335
336#endif // IFPACK2_DETAILS_MULTIVECTORLOCALGATHERSCATTER_HPP
Implementation detail of Ifpack2::Container subclasses.
Definition Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74