Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
45 #include "Stokhos_Sacado_Kokkos_MP_Vector.hpp"
46 #include "Tpetra_Map.hpp"
47 #include "Tpetra_MultiVector.hpp"
48 #include "Tpetra_CrsGraph.hpp"
49 #include "Tpetra_CrsMatrix.hpp"
50 
51 namespace 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> >
97  create_flat_mp_graph(
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> > >
171  create_flat_mp_graph(
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,
250  LocalOrdinal,GlobalOrdinal,Node> >
251  create_flat_vector_view(
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 
259  typedef Sacado::MP::Vector<Storage> Scalar;
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,
297  LocalOrdinal,GlobalOrdinal,Node> >
298  create_flat_vector_view(
299  const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
300  LocalOrdinal,GlobalOrdinal,Node>& vec,
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,
315  LocalOrdinal,GlobalOrdinal,Node> >
316  create_flat_vector_view(
317  Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
318  LocalOrdinal,GlobalOrdinal,Node>& vec,
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 
324  typedef Sacado::MP::Vector<Storage> Scalar;
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,
357  LocalOrdinal,GlobalOrdinal,Node> >
358  create_flat_vector_view(
359  Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
360  LocalOrdinal,GlobalOrdinal,Node>& vec,
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,
375  LocalOrdinal,GlobalOrdinal,Node> >
376  create_flat_vector_view(
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,
391  LocalOrdinal,GlobalOrdinal,Node> >
392  create_flat_vector_view(
393  const Tpetra::Vector<Sacado::MP::Vector<Storage>,
394  LocalOrdinal,GlobalOrdinal,Node>& vec,
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,
409  LocalOrdinal,GlobalOrdinal,Node> >
410  create_flat_vector_view(
411  Tpetra::Vector<Sacado::MP::Vector<Storage>,
412  LocalOrdinal,GlobalOrdinal,Node>& vec,
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,
425  LocalOrdinal,GlobalOrdinal,Node> >
426  create_flat_vector_view(
427  Tpetra::Vector<Sacado::MP::Vector<Storage>,
428  LocalOrdinal,GlobalOrdinal,Node>& vec,
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,
444  LocalOrdinal,GlobalOrdinal,
445  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
446  create_flat_vector_view(
447  const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
448  LocalOrdinal,GlobalOrdinal,
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,
481  LocalOrdinal,GlobalOrdinal,
482  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
483  create_flat_vector_view(
484  Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
485  LocalOrdinal,GlobalOrdinal,
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,
512  LocalOrdinal,GlobalOrdinal,Node> >
513  create_flat_matrix(
514  const Tpetra::CrsMatrix<Sacado::MP::Vector<Storage>,
515  LocalOrdinal,GlobalOrdinal,Node>& mat,
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 
523  typedef Sacado::MP::Vector<Storage> Scalar;
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
Storage
Stokhos::StandardStorage< int, double > Storage
Definition: linear2d_diffusion_scalar_types.hpp:47
Stokhos_Sacado_Kokkos_MP_Vector.hpp
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Scalar
pce_type Scalar
Definition: linear2d_diffusion_scalar_types.hpp:49
Stokhos::create_flat_matrix
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)
Definition: Stokhos_Tpetra_Utilities_UQ_PCE.hpp:439
LocalOrdinal
int LocalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:50
Stokhos
Top-level namespace for Stokhos classes and functions.
Definition: Stokhos_AbstractPreconditionerFactory.hpp:48
Node
KokkosClassic::DefaultNode::DefaultNodeType Node
Definition: experimental/linear2d_diffusion_pce.cpp:79
Kokkos::deep_copy
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
Definition: Stokhos_CrsMatrix.hpp:688
Stokhos::create_flat_mp_graph
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)
Definition: Stokhos_Tpetra_Utilities_MP_Vector.hpp:97
GlobalOrdinal
int GlobalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:51
Stokhos::create_flat_vector_view
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)
Definition: Stokhos_Tpetra_Utilities_UQ_PCE.hpp:238
Stokhos::create_flat_map
Teuchos::RCP< Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > create_flat_map(const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map, const LocalOrdinal block_size)
Definition: Stokhos_Tpetra_Utilities_MP_Vector.hpp:57
Kokkos::create_mirror_view
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
Definition: Stokhos_CrsMatrix.hpp:677
Generated by   doxygen 1.8.14