Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_Tpetra_Utilities_MP_Vector.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos 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 Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP
43#define STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP
44
46#include "Tpetra_Map.hpp"
47#include "Tpetra_MultiVector.hpp"
48#include "Tpetra_CrsGraph.hpp"
49#include "Tpetra_CrsMatrix.hpp"
50
51namespace Stokhos {
52
53 // Create a flattened map for a map representing a distribution for an
54 // embedded scalar type
55 template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
56 Teuchos::RCP< Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
57 create_flat_map(const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map,
58 const LocalOrdinal block_size) {
59 using Tpetra::global_size_t;
60 using Teuchos::ArrayView;
61 using Teuchos::Array;
62 using Teuchos::RCP;
63 using Teuchos::rcp;
64
65 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
66
67 // Get map info
68 const global_size_t num_global_entries = map.getGlobalNumElements();
69 const size_t num_local_entries = map.getLocalNumElements();
70 const GlobalOrdinal index_base = map.getIndexBase();
71 ArrayView<const GlobalOrdinal> element_list =
72 map.getLocalElementList();
73
74 // Create new elements
75 const global_size_t flat_num_global_entries = num_global_entries*block_size;
76 const size_t flat_num_local_entries = num_local_entries * block_size;
77 const GlobalOrdinal flat_index_base = index_base;
78 Array<GlobalOrdinal> flat_element_list(flat_num_local_entries);
79 for (size_t i=0; i<num_local_entries; ++i)
80 for (LocalOrdinal j=0; j<block_size; ++j)
81 flat_element_list[i*block_size+j] = element_list[i]*block_size+j;
82
83 // Create new map
84 RCP<Map> flat_map =
85 rcp(new Map(flat_num_global_entries, flat_element_list(),
86 flat_index_base, map.getComm()));
87
88 return flat_map;
89 }
90
91 // Create a flattened graph for a graph from a matrix with the
92 // MP::Vector scalar type (each block is an identity matrix)
93 // If flat_domain_map and/or flat_range_map are null, they will be computed,
94 // otherwise they will be used as-is.
95 template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
96 Teuchos::RCP< Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >
98 const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>& graph,
99 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_domain_map,
100 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_range_map,
101 const LocalOrdinal block_size) {
102 using Teuchos::ArrayRCP;
103 using Teuchos::RCP;
104 using Teuchos::rcp;
105
106 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
107 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
108
109 // Build domain map if necessary
110 if (flat_domain_map == Teuchos::null)
111 flat_domain_map = create_flat_map(*(graph.getDomainMap()), block_size);
112
113 // Build range map if necessary
114 if (flat_range_map == Teuchos::null)
115 flat_range_map = create_flat_map(*(graph.getRangeMap()), block_size);
116
117 // Build column map
118 RCP<const Map> flat_col_map =
119 create_flat_map(*(graph.getColMap()), block_size);
120
121 // Build row map if necessary
122 // Check if range_map == row_map, then we can use flat_range_map
123 // as the flattened row map
124 RCP<const Map> flat_row_map;
125 if (graph.getRangeMap() == graph.getRowMap())
126 flat_row_map = flat_range_map;
127 else
128 flat_row_map = create_flat_map(*(graph.getRowMap()), block_size);
129
130 // Build flattened row offsets and column indices
131 auto row_offsets = graph.getLocalRowPtrsHost();
132 auto col_indices = graph.getLocalIndicesHost();
133 const size_t num_row = graph.getLocalNumRows();
134 const size_t num_col_indices = col_indices.size();
135 ArrayRCP<size_t> flat_row_offsets(num_row*block_size+1);
136 ArrayRCP<LocalOrdinal> flat_col_indices(num_col_indices * block_size);
137 for (size_t row=0; row<num_row; ++row) {
138 const size_t row_beg = row_offsets[row];
139 const size_t row_end = row_offsets[row+1];
140 const size_t num_col = row_end - row_beg;
141 for (LocalOrdinal j=0; j<block_size; ++j) {
142 const size_t flat_row = row*block_size + j;
143 const size_t flat_row_beg = row_beg*block_size + j*num_col;
144 flat_row_offsets[flat_row] = flat_row_beg;
145 for (size_t entry=0; entry<num_col; ++entry) {
146 const LocalOrdinal col = col_indices[row_beg+entry];
147 const LocalOrdinal flat_col = col*block_size + j;
148 flat_col_indices[flat_row_beg+entry] = flat_col;
149 }
150 }
151 }
152 flat_row_offsets[num_row*block_size] = num_col_indices*block_size;
153
154 // Build flattened graph
155 RCP<Graph> flat_graph =
156 rcp(new Graph(flat_row_map, flat_col_map,
157 flat_row_offsets, flat_col_indices));
158 flat_graph->fillComplete(flat_domain_map, flat_range_map);
159
160 return flat_graph;
161 }
162
163
164 // Create a flattened graph for a graph from a matrix with the
165 // MP::Vector scalar type (each block is an identity matrix)
166 // If flat_domain_map and/or flat_range_map are null, they will be computed,
167 // otherwise they will be used as-is.
168 template <typename LocalOrdinal, typename GlobalOrdinal, typename Device>
169 Teuchos::RCP< Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,
170 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
172 const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& graph,
173 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_domain_map,
174 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_range_map,
175 const LocalOrdinal block_size) {
176 using Teuchos::ArrayRCP;
177 using Teuchos::RCP;
178 using Teuchos::rcp;
179
180 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
181 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
182 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
183 typedef typename Graph::local_graph_device_type::row_map_type::non_const_type RowPtrs;
184 typedef typename Graph::local_graph_device_type::entries_type::non_const_type LocalIndices;
185
186 // Build domain map if necessary
187 if (flat_domain_map == Teuchos::null)
188 flat_domain_map = create_flat_map(*(graph.getDomainMap()), block_size);
189
190 // Build range map if necessary
191 if (flat_range_map == Teuchos::null)
192 flat_range_map = create_flat_map(*(graph.getRangeMap()), block_size);
193
194 // Build column map
195 RCP<const Map> flat_col_map =
196 create_flat_map(*(graph.getColMap()), block_size);
197
198 // Build row map if necessary
199 // Check if range_map == row_map, then we can use flat_range_map
200 // as the flattened row map
201 RCP<const Map> flat_row_map;
202 if (graph.getRangeMap() == graph.getRowMap())
203 flat_row_map = flat_range_map;
204 else
205 flat_row_map = create_flat_map(*(graph.getRowMap()), block_size);
206
207 // Build flattened row offsets and column indices
208 auto row_offsets = graph.getLocalRowPtrsHost();
209 auto col_indices = graph.getLocalIndicesHost();
210 const size_t num_row = graph.getLocalNumRows();
211 const size_t num_col_indices = col_indices.size();
212 RowPtrs flat_row_offsets("row_ptrs", num_row*block_size+1);
213 LocalIndices flat_col_indices("col_indices", num_col_indices * block_size);
214 auto flat_row_offsets_host = Kokkos::create_mirror_view(flat_row_offsets);
215 auto flat_col_indices_host = Kokkos::create_mirror_view(flat_col_indices);
216 for (size_t row=0; row<num_row; ++row) {
217 const size_t row_beg = row_offsets[row];
218 const size_t row_end = row_offsets[row+1];
219 const size_t num_col = row_end - row_beg;
220 for (LocalOrdinal j=0; j<block_size; ++j) {
221 const size_t flat_row = row*block_size + j;
222 const size_t flat_row_beg = row_beg*block_size + j*num_col;
223 flat_row_offsets_host[flat_row] = flat_row_beg;
224 for (size_t entry=0; entry<num_col; ++entry) {
225 const LocalOrdinal col = col_indices[row_beg+entry];
226 const LocalOrdinal flat_col = col*block_size + j;
227 flat_col_indices_host[flat_row_beg+entry] = flat_col;
228 }
229 }
230 }
231 flat_row_offsets_host[num_row*block_size] = num_col_indices*block_size;
232 Kokkos::deep_copy(flat_row_offsets, flat_row_offsets_host);
233 Kokkos::deep_copy(flat_col_indices, flat_col_indices_host);
234
235 // Build flattened graph
236 RCP<Graph> flat_graph =
237 rcp(new Graph(flat_row_map, flat_col_map,
238 flat_row_offsets, flat_col_indices));
239 flat_graph->fillComplete(flat_domain_map, flat_range_map);
240
241 return flat_graph;
242 }
243
244
245 // Create a flattened vector by unrolling the MP::Vector scalar type. The
246 // returned vector is a view of the original
247 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
248 typename Node>
249 Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
252 const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
253 LocalOrdinal,GlobalOrdinal,Node>& vec_const,
254 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
255 using Teuchos::ArrayRCP;
256 using Teuchos::RCP;
257 using Teuchos::rcp;
258
260 typedef typename Storage::value_type BaseScalar;
261 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Vector;
262 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
263
264 // MP size
265 const LocalOrdinal mp_size = Storage::static_size;
266
267 // Cast-away const since resulting vector is a view and we can't create
268 // a const-vector directly
269 Vector& vec = const_cast<Vector&>(vec_const);
270
271 // Get data
272 ArrayRCP<Scalar> vec_vals = vec.get1dViewNonConst();
273 const size_t vec_size = vec_vals.size();
274
275 // Create view of data
276 BaseScalar *flat_vec_ptr =
277 reinterpret_cast<BaseScalar*>(vec_vals.getRawPtr());
278 ArrayRCP<BaseScalar> flat_vec_vals =
279 Teuchos::arcp(flat_vec_ptr, 0, vec_size*mp_size, false);
280
281 // Create flat vector
282 const size_t stride = vec.getStride();
283 const size_t flat_stride = stride * mp_size;
284 const size_t num_vecs = vec.getNumVectors();
285 RCP<FlatVector> flat_vec =
286 rcp(new FlatVector(flat_map, flat_vec_vals, flat_stride, num_vecs));
287
288 return flat_vec;
289 }
290
291 // Create a flattened vector by unrolling the MP::Vector scalar type. The
292 // returned vector is a view of the original. This version creates the
293 // map if necessary
294 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
295 typename Node>
296 Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
299 const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
301 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
302 if (flat_map == Teuchos::null) {
303 const LocalOrdinal mp_size = Storage::static_size;
304 flat_map = create_flat_map(*(vec.getMap()), mp_size);
305 }
306 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
307 return create_flat_vector_view(vec, const_flat_map);
308 }
309
310 // Create a flattened vector by unrolling the MP::Vector scalar type. The
311 // returned vector is a view of the original
312 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
313 typename Node>
314 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
317 Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
319 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
320 using Teuchos::ArrayRCP;
321 using Teuchos::RCP;
322 using Teuchos::rcp;
323
325 typedef typename Storage::value_type BaseScalar;
326 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
327
328 // MP size
329 const LocalOrdinal mp_size = Storage::static_size;
330
331 // Get data
332 ArrayRCP<Scalar> vec_vals = vec.get1dViewNonConst();
333 const size_t vec_size = vec_vals.size();
334
335 // Create view of data
336 BaseScalar *flat_vec_ptr =
337 reinterpret_cast<BaseScalar*>(vec_vals.getRawPtr());
338 ArrayRCP<BaseScalar> flat_vec_vals =
339 Teuchos::arcp(flat_vec_ptr, 0, vec_size*mp_size, false);
340
341 // Create flat vector
342 const size_t stride = vec.getStride();
343 const size_t flat_stride = stride * mp_size;
344 const size_t num_vecs = vec.getNumVectors();
345 RCP<FlatVector> flat_vec =
346 rcp(new FlatVector(flat_map, flat_vec_vals, flat_stride, num_vecs));
347
348 return flat_vec;
349 }
350
351 // Create a flattened vector by unrolling the MP::Vector scalar type. The
352 // returned vector is a view of the original. This version creates the
353 // map if necessary
354 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
355 typename Node>
356 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
359 Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
361 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
362 if (flat_map == Teuchos::null) {
363 const LocalOrdinal mp_size = Storage::static_size;
364 flat_map = create_flat_map(*(vec.getMap()), mp_size);
365 }
366 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
367 return create_flat_vector_view(vec, const_flat_map);
368 }
369
370 // Create a flattened vector by unrolling the MP::Vector scalar type. The
371 // returned vector is a view of the original
372 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
373 typename Node>
374 Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
377 const Tpetra::Vector<Sacado::MP::Vector<Storage>,
378 LocalOrdinal,GlobalOrdinal,Node>& vec_const,
379 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
380 const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec_const;
381 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv = create_flat_vector_view(mv, flat_map);
382 return flat_mv->getVector(0);
383 }
384
385 // Create a flattened vector by unrolling the MP::Vector scalar type. The
386 // returned vector is a view of the original. This version creates the
387 // map if necessary
388 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
389 typename Node>
390 Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
393 const Tpetra::Vector<Sacado::MP::Vector<Storage>,
395 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
396 if (flat_map == Teuchos::null) {
397 const LocalOrdinal mp_size = Storage::static_size;
398 flat_map = create_flat_map(*(vec.getMap()), mp_size);
399 }
400 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
401 return create_flat_vector_view(vec, const_flat_map);
402 }
403
404 // Create a flattened vector by unrolling the MP::Vector scalar type. The
405 // returned vector is a view of the original
406 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
407 typename Node>
408 Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
411 Tpetra::Vector<Sacado::MP::Vector<Storage>,
413 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
414 Tpetra::MultiVector<Sacado::MP::Vector<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec;
415 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv = create_flat_vector_view(mv, flat_map);
416 return flat_mv->getVectorNonConst(0);
417 }
418
419 // Create a flattened vector by unrolling the MP::Vector scalar type. The
420 // returned vector is a view of the original. This version creates the
421 // map if necessary
422 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
423 typename Node>
424 Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
427 Tpetra::Vector<Sacado::MP::Vector<Storage>,
429 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
430 if (flat_map == Teuchos::null) {
431 const LocalOrdinal mp_size = Storage::static_size;
432 flat_map = create_flat_map(*(vec.getMap()), mp_size);
433 }
434 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
435 return create_flat_vector_view(vec, const_flat_map);
436 }
437
438
439 // Create a flattened vector by unrolling the MP::Vector scalar type. The
440 // returned vector is a view of the original
441 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
442 typename Device>
443 Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
445 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
447 const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
449 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
450 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
451 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
452 using Teuchos::RCP;
453 using Teuchos::rcp;
454
455 typedef typename Storage::value_type BaseScalar;
456 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
457 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
458 typedef typename FlatVector::dual_view_type::t_dev flat_view_type;
459
460 // Have to do a nasty const-cast because getLocalViewDevice(ReadWrite) is a
461 // non-const method, yet getLocalViewDevice(ReadOnly) returns a const-view
462 // (i.e., with a constant scalar type), and there is no way to make a
463 // MultiVector out of it!
464 typedef Tpetra::MultiVector<Sacado::MP::Vector<Storage>, LocalOrdinal,GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode<Device> > mv_type;
465 mv_type& vec_nc = const_cast<mv_type&>(vec);
466
467 // Create flattenend view using special reshaping view assignment operator
468 flat_view_type flat_vals = vec_nc.getLocalViewDevice(Tpetra::Access::ReadWrite);
469
470 // Create flat vector
471 RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
472
473 return flat_vec;
474 }
475
476 // Create a flattened vector by unrolling the MP::Vector scalar type. The
477 // returned vector is a view of the original
478 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
479 typename Device>
480 Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
482 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
484 Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
486 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
487 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
488 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
489 using Teuchos::RCP;
490 using Teuchos::rcp;
491
492 typedef typename Storage::value_type BaseScalar;
493 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
494 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
495 typedef typename FlatVector::dual_view_type::t_dev flat_view_type;
496
497 // Create flattenend view using special reshaping view assignment operator
498 flat_view_type flat_vals = vec.getLocalViewDevice(Tpetra::Access::ReadWrite);
499
500 // Create flat vector
501 RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
502
503 return flat_vec;
504 }
505
506
507 // Create a flattened matrix by unrolling the MP::Vector scalar type. The
508 // returned matrix is NOT a view of the original (and can't be)
509 template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
510 typename Node>
511 Teuchos::RCP< Tpetra::CrsMatrix<typename Storage::value_type,
514 const Tpetra::CrsMatrix<Sacado::MP::Vector<Storage>,
516 const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& flat_graph,
517 const LocalOrdinal block_size) {
518 using Teuchos::ArrayView;
519 using Teuchos::Array;
520 using Teuchos::RCP;
521 using Teuchos::rcp;
522
524 typedef typename Storage::value_type BaseScalar;
525 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Matrix;
526 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
527
528 // Create flat matrix
529 RCP<FlatMatrix> flat_mat = rcp(new FlatMatrix(flat_graph));
530
531 // Set values
532 const size_t num_rows = mat.getLocalNumRows();
533 const size_t max_cols = mat.getLocalMaxNumRowEntries();
534 typename Matrix::local_inds_host_view_type indices, flat_indices;
535 typename Matrix::values_host_view_type values;
536 Array<BaseScalar> flat_values(max_cols);
537 for (size_t row=0; row<num_rows; ++row) {
538 mat.getLocalRowView(row, indices, values);
539 const size_t num_col = mat.getNumEntriesInLocalRow(row);
540 for (LocalOrdinal i=0; i<block_size; ++i) {
541 const LocalOrdinal flat_row = row*block_size + i;
542 for (size_t j=0; j<num_col; ++j)
543 flat_values[j] = values[j].coeff(i);
544 flat_graph->getLocalRowView(flat_row, flat_indices);
545 flat_mat->replaceLocalValues(flat_row, Kokkos::Compat::getConstArrayView(flat_indices),
546 flat_values(0, num_col));
547 }
548 }
549 flat_mat->fillComplete(flat_graph->getDomainMap(),
550 flat_graph->getRangeMap());
551
552 return flat_mat;
553 }
554
555} // namespace Stokhos
556
557#endif // STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP
KokkosClassic::DefaultNode::DefaultNodeType Node
Stokhos::StandardStorage< int, double > Storage
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
Top-level namespace for Stokhos classes and functions.
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_map)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &cijk_graph, const CijkType &cijk_dev)
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > create_flat_mp_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_range_map, const LocalOrdinal block_size)
Teuchos::RCP< Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > create_flat_map(const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map, const LocalOrdinal block_size)