Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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.getNodeNumElements();
70  const GlobalOrdinal index_base = map.getIndexBase();
71  ArrayView<const GlobalOrdinal> element_list =
72  map.getNodeElementList();
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  ArrayRCP<const size_t> row_offsets = graph.getNodeRowPtrs();
132  ArrayRCP<const LocalOrdinal> col_indices = graph.getNodePackedIndices();
133  const size_t num_row = graph.getNodeNumRows();
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 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
164 
165  // Create a flattened graph for a graph from a matrix with the
166  // MP::Vector scalar type (each block is an identity matrix)
167  // If flat_domain_map and/or flat_range_map are null, they will be computed,
168  // otherwise they will be used as-is.
169  template <typename LocalOrdinal, typename GlobalOrdinal, typename Device>
170  Teuchos::RCP< Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,
171  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
172  create_flat_mp_graph(
173  const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& graph,
174  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_domain_map,
175  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_range_map,
176  const LocalOrdinal block_size) {
177  using Teuchos::ArrayRCP;
178  using Teuchos::RCP;
179  using Teuchos::rcp;
180 
181  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
182  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
183  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
184  typedef typename Graph::local_graph_type::row_map_type::non_const_type RowPtrs;
185  typedef typename Graph::local_graph_type::entries_type::non_const_type LocalIndices;
186 
187  // Build domain map if necessary
188  if (flat_domain_map == Teuchos::null)
189  flat_domain_map = create_flat_map(*(graph.getDomainMap()), block_size);
190 
191  // Build range map if necessary
192  if (flat_range_map == Teuchos::null)
193  flat_range_map = create_flat_map(*(graph.getRangeMap()), block_size);
194 
195  // Build column map
196  RCP<const Map> flat_col_map =
197  create_flat_map(*(graph.getColMap()), block_size);
198 
199  // Build row map if necessary
200  // Check if range_map == row_map, then we can use flat_range_map
201  // as the flattened row map
202  RCP<const Map> flat_row_map;
203  if (graph.getRangeMap() == graph.getRowMap())
204  flat_row_map = flat_range_map;
205  else
206  flat_row_map = create_flat_map(*(graph.getRowMap()), block_size);
207 
208  // Build flattened row offsets and column indices
209  ArrayRCP<const size_t> row_offsets = graph.getNodeRowPtrs();
210  ArrayRCP<const LocalOrdinal> col_indices = graph.getNodePackedIndices();
211  const size_t num_row = graph.getNodeNumRows();
212  const size_t num_col_indices = col_indices.size();
213  RowPtrs flat_row_offsets("row_ptrs", num_row*block_size+1);
214  LocalIndices flat_col_indices("col_indices", num_col_indices * block_size);
215  for (size_t row=0; row<num_row; ++row) {
216  const size_t row_beg = row_offsets[row];
217  const size_t row_end = row_offsets[row+1];
218  const size_t num_col = row_end - row_beg;
219  for (LocalOrdinal j=0; j<block_size; ++j) {
220  const size_t flat_row = row*block_size + j;
221  const size_t flat_row_beg = row_beg*block_size + j*num_col;
222  flat_row_offsets[flat_row] = flat_row_beg;
223  for (size_t entry=0; entry<num_col; ++entry) {
224  const LocalOrdinal col = col_indices[row_beg+entry];
225  const LocalOrdinal flat_col = col*block_size + j;
226  flat_col_indices[flat_row_beg+entry] = flat_col;
227  }
228  }
229  }
230  flat_row_offsets[num_row*block_size] = num_col_indices*block_size;
231 
232  // Build flattened graph
233  RCP<Graph> flat_graph =
234  rcp(new Graph(flat_row_map, flat_col_map,
235  flat_row_offsets, flat_col_indices));
236  flat_graph->fillComplete(flat_domain_map, flat_range_map);
237 
238  return flat_graph;
239  }
240 
241 #endif
242 
243  // Create a flattened vector by unrolling the MP::Vector scalar type. The
244  // returned vector is a view of the original
245  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
246  typename Node>
247  Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
248  LocalOrdinal,GlobalOrdinal,Node> >
249  create_flat_vector_view(
250  const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
251  LocalOrdinal,GlobalOrdinal,Node>& vec_const,
252  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
253  using Teuchos::ArrayRCP;
254  using Teuchos::RCP;
255  using Teuchos::rcp;
256 
257  typedef Sacado::MP::Vector<Storage> Scalar;
258  typedef typename Storage::value_type BaseScalar;
259  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Vector;
260  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
261 
262  // MP size
263  const LocalOrdinal mp_size = Storage::static_size;
264 
265  // Cast-away const since resulting vector is a view and we can't create
266  // a const-vector directly
267  Vector& vec = const_cast<Vector&>(vec_const);
268 
269  // Get data
270  ArrayRCP<Scalar> vec_vals = vec.get1dViewNonConst();
271  const size_t vec_size = vec_vals.size();
272 
273  // Create view of data
274  BaseScalar *flat_vec_ptr =
275  reinterpret_cast<BaseScalar*>(vec_vals.getRawPtr());
276  ArrayRCP<BaseScalar> flat_vec_vals =
277  Teuchos::arcp(flat_vec_ptr, 0, vec_size*mp_size, false);
278 
279  // Create flat vector
280  const size_t stride = vec.getStride();
281  const size_t flat_stride = stride * mp_size;
282  const size_t num_vecs = vec.getNumVectors();
283  RCP<FlatVector> flat_vec =
284  rcp(new FlatVector(flat_map, flat_vec_vals, flat_stride, num_vecs));
285 
286  return flat_vec;
287  }
288 
289  // Create a flattened vector by unrolling the MP::Vector scalar type. The
290  // returned vector is a view of the original. This version creates the
291  // map if necessary
292  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
293  typename Node>
294  Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
295  LocalOrdinal,GlobalOrdinal,Node> >
296  create_flat_vector_view(
297  const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
298  LocalOrdinal,GlobalOrdinal,Node>& vec,
299  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
300  if (flat_map == Teuchos::null) {
301  const LocalOrdinal mp_size = Storage::static_size;
302  flat_map = create_flat_map(*(vec.getMap()), mp_size);
303  }
304  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
305  return create_flat_vector_view(vec, const_flat_map);
306  }
307 
308  // Create a flattened vector by unrolling the MP::Vector scalar type. The
309  // returned vector is a view of the original
310  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
311  typename Node>
312  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
313  LocalOrdinal,GlobalOrdinal,Node> >
314  create_flat_vector_view(
315  Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
316  LocalOrdinal,GlobalOrdinal,Node>& vec,
317  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
318  using Teuchos::ArrayRCP;
319  using Teuchos::RCP;
320  using Teuchos::rcp;
321 
322  typedef Sacado::MP::Vector<Storage> Scalar;
323  typedef typename Storage::value_type BaseScalar;
324  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
325 
326  // MP size
327  const LocalOrdinal mp_size = Storage::static_size;
328 
329  // Get data
330  ArrayRCP<Scalar> vec_vals = vec.get1dViewNonConst();
331  const size_t vec_size = vec_vals.size();
332 
333  // Create view of data
334  BaseScalar *flat_vec_ptr =
335  reinterpret_cast<BaseScalar*>(vec_vals.getRawPtr());
336  ArrayRCP<BaseScalar> flat_vec_vals =
337  Teuchos::arcp(flat_vec_ptr, 0, vec_size*mp_size, false);
338 
339  // Create flat vector
340  const size_t stride = vec.getStride();
341  const size_t flat_stride = stride * mp_size;
342  const size_t num_vecs = vec.getNumVectors();
343  RCP<FlatVector> flat_vec =
344  rcp(new FlatVector(flat_map, flat_vec_vals, flat_stride, num_vecs));
345 
346  return flat_vec;
347  }
348 
349  // Create a flattened vector by unrolling the MP::Vector scalar type. The
350  // returned vector is a view of the original. This version creates the
351  // map if necessary
352  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
353  typename Node>
354  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
355  LocalOrdinal,GlobalOrdinal,Node> >
356  create_flat_vector_view(
357  Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
358  LocalOrdinal,GlobalOrdinal,Node>& vec,
359  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
360  if (flat_map == Teuchos::null) {
361  const LocalOrdinal mp_size = Storage::static_size;
362  flat_map = create_flat_map(*(vec.getMap()), mp_size);
363  }
364  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
365  return create_flat_vector_view(vec, const_flat_map);
366  }
367 
368  // Create a flattened vector by unrolling the MP::Vector scalar type. The
369  // returned vector is a view of the original
370  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
371  typename Node>
372  Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
373  LocalOrdinal,GlobalOrdinal,Node> >
374  create_flat_vector_view(
375  const Tpetra::Vector<Sacado::MP::Vector<Storage>,
376  LocalOrdinal,GlobalOrdinal,Node>& vec_const,
377  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
378  const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec_const;
379  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv = create_flat_vector_view(mv, flat_map);
380  return flat_mv->getVector(0);
381  }
382 
383  // Create a flattened vector by unrolling the MP::Vector scalar type. The
384  // returned vector is a view of the original. This version creates the
385  // map if necessary
386  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
387  typename Node>
388  Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
389  LocalOrdinal,GlobalOrdinal,Node> >
390  create_flat_vector_view(
391  const Tpetra::Vector<Sacado::MP::Vector<Storage>,
392  LocalOrdinal,GlobalOrdinal,Node>& vec,
393  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
394  if (flat_map == Teuchos::null) {
395  const LocalOrdinal mp_size = Storage::static_size;
396  flat_map = create_flat_map(*(vec.getMap()), mp_size);
397  }
398  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
399  return create_flat_vector_view(vec, const_flat_map);
400  }
401 
402  // Create a flattened vector by unrolling the MP::Vector scalar type. The
403  // returned vector is a view of the original
404  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
405  typename Node>
406  Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
407  LocalOrdinal,GlobalOrdinal,Node> >
408  create_flat_vector_view(
409  Tpetra::Vector<Sacado::MP::Vector<Storage>,
410  LocalOrdinal,GlobalOrdinal,Node>& vec,
411  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
412  Tpetra::MultiVector<Sacado::MP::Vector<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec;
413  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv = create_flat_vector_view(mv, flat_map);
414  return flat_mv->getVectorNonConst(0);
415  }
416 
417  // Create a flattened vector by unrolling the MP::Vector scalar type. The
418  // returned vector is a view of the original. This version creates the
419  // map if necessary
420  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
421  typename Node>
422  Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
423  LocalOrdinal,GlobalOrdinal,Node> >
424  create_flat_vector_view(
425  Tpetra::Vector<Sacado::MP::Vector<Storage>,
426  LocalOrdinal,GlobalOrdinal,Node>& vec,
427  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& flat_map) {
428  if (flat_map == Teuchos::null) {
429  const LocalOrdinal mp_size = Storage::static_size;
430  flat_map = create_flat_map(*(vec.getMap()), mp_size);
431  }
432  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
433  return create_flat_vector_view(vec, const_flat_map);
434  }
435 
436 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
437 
438  // Create a flattened vector by unrolling the MP::Vector scalar type. The
439  // returned vector is a view of the original
440  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
441  typename Device>
442  Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
443  LocalOrdinal,GlobalOrdinal,
444  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
445  create_flat_vector_view(
446  const Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
447  LocalOrdinal,GlobalOrdinal,
448  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
449  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
450  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
451  using Teuchos::RCP;
452  using Teuchos::rcp;
453 
454  typedef typename Storage::value_type BaseScalar;
455  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
456  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
457  typedef typename FlatVector::dual_view_type flat_view_type;
458 
459  // Create flattenend view using special reshaping view assignment operator
460  flat_view_type flat_vals (vec.getLocalViewDevice(), vec.getLocalViewHost());
461  if (vec.need_sync_device ()) {
462  flat_vals.modify_host ();
463  }
464  else if (vec.need_sync_host ()) {
465  flat_vals.modify_device ();
466  }
467 
468  // Create flat vector
469  RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
470 
471  return flat_vec;
472  }
473 
474  // Create a flattened vector by unrolling the MP::Vector scalar type. The
475  // returned vector is a view of the original
476  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
477  typename Device>
478  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
479  LocalOrdinal,GlobalOrdinal,
480  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
481  create_flat_vector_view(
482  Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
483  LocalOrdinal,GlobalOrdinal,
484  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
485  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
486  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
487  using Teuchos::RCP;
488  using Teuchos::rcp;
489 
490  typedef typename Storage::value_type BaseScalar;
491  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
492  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
493  typedef typename FlatVector::dual_view_type flat_view_type;
494 
495  // Create flattenend view using special reshaping view assignment operator
496  flat_view_type flat_vals (vec.getLocalViewDevice(), vec.getLocalViewHost());
497  if (vec.need_sync_device ()) {
498  flat_vals.modify_host ();
499  }
500  else if (vec.need_sync_host ()) {
501  flat_vals.modify_device ();
502  }
503 
504  // Create flat vector
505  RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
506 
507  return flat_vec;
508  }
509 
510 #endif
511 
512  // Create a flattened matrix by unrolling the MP::Vector scalar type. The
513  // returned matrix is NOT a view of the original (and can't be)
514  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
515  typename Node>
516  Teuchos::RCP< Tpetra::CrsMatrix<typename Storage::value_type,
517  LocalOrdinal,GlobalOrdinal,Node> >
518  create_flat_matrix(
519  const Tpetra::CrsMatrix<Sacado::MP::Vector<Storage>,
520  LocalOrdinal,GlobalOrdinal,Node>& mat,
521  const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& flat_graph,
522  const LocalOrdinal block_size) {
523  using Teuchos::ArrayView;
524  using Teuchos::Array;
525  using Teuchos::RCP;
526  using Teuchos::rcp;
527 
528  typedef Sacado::MP::Vector<Storage> Scalar;
529  typedef typename Storage::value_type BaseScalar;
530  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
531 
532  // Create flat matrix
533  RCP<FlatMatrix> flat_mat = rcp(new FlatMatrix(flat_graph));
534 
535  // Set values
536  const size_t num_rows = mat.getNodeNumRows();
537  const size_t max_cols = mat.getNodeMaxNumRowEntries();
538  ArrayView<const LocalOrdinal> indices, flat_indices;
539  ArrayView<const Scalar> values;
540  Array<BaseScalar> flat_values(max_cols);
541  for (size_t row=0; row<num_rows; ++row) {
542  mat.getLocalRowView(row, indices, values);
543  const size_t num_col = mat.getNumEntriesInLocalRow(row);
544  for (LocalOrdinal i=0; i<block_size; ++i) {
545  const LocalOrdinal flat_row = row*block_size + i;
546  for (size_t j=0; j<num_col; ++j)
547  flat_values[j] = values[j].coeff(i);
548  flat_graph->getLocalRowView(flat_row, flat_indices);
549  flat_mat->replaceLocalValues(flat_row, flat_indices,
550  flat_values(0, num_col));
551  }
552  }
553  flat_mat->fillComplete(flat_graph->getDomainMap(),
554  flat_graph->getRangeMap());
555 
556  return flat_mat;
557  }
558 
559 } // namespace Stokhos
560 
561 #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
Teuchos::ArrayRCP::getRawPtr
T * getRawPtr() const
Teuchos::ArrayRCP::size
size_type size() const
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Scalar
pce_type Scalar
Definition: linear2d_diffusion_scalar_types.hpp:49
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Stokhos::create_flat_vector_view
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec_const, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
Definition: Stokhos_Tpetra_Utilities_MP_Vector.hpp:249
LocalOrdinal
int LocalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:50
Node
KokkosClassic::DefaultNode::DefaultNodeType Node
Definition: experimental/linear2d_diffusion_pce.cpp:79
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
Stokhos::create_flat_matrix
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &flat_graph, const LocalOrdinal block_size)
Definition: Stokhos_Tpetra_Utilities_MP_Vector.hpp:518
GlobalOrdinal
int GlobalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:51
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
Generated by   doxygen 1.8.5