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_UQ_PCE.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_UQ_PCE_HPP
43 #define STOKHOS_TPETRA_UTILITIES_UQ_PCE_HPP
44 
45 #include "Stokhos_Sacado_Kokkos_UQ_PCE.hpp"
46 #include "Stokhos_Tpetra_Utilities_MP_Vector.hpp"
47 #include "Tpetra_Map.hpp"
48 #include "Tpetra_MultiVector.hpp"
49 #include "Tpetra_CrsGraph.hpp"
50 #include "Tpetra_CrsMatrix.hpp"
51 
52 namespace Stokhos {
53 
54 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
55 
56  // Build a CRS graph from a sparse Cijk tensor
57  template <typename LocalOrdinal, typename GlobalOrdinal, typename Device,
58  typename CijkType>
59  Teuchos::RCP< Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,
60  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
61  create_cijk_crs_graph(const CijkType& cijk,
62  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
63  const size_t matrix_pce_size) {
64  using Teuchos::RCP;
65  using Teuchos::arrayView;
66 
67  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
68  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
69  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
70 
71  const size_t pce_sz = cijk.dimension();
72  RCP<const Map> map =
73  Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(pce_sz, comm);
74  RCP<Graph> graph;
75  if (matrix_pce_size == 1) {
76  graph = Tpetra::createCrsGraph(map, 1);
77  // Mean-based case -- graph is diagonal
78  for (size_t i=0; i<pce_sz; ++i) {
79  const GlobalOrdinal row = i;
80  graph->insertGlobalIndices(row, arrayView(&row, 1));
81  }
82  }
83  else {
84  // General case
85 
86  // Get max num entries
87  size_t max_num_entry = 0;
88  for (size_t i=0; i<pce_sz; ++i) {
89  const size_t num_entry = cijk.num_entry(i);
90  max_num_entry = (num_entry > max_num_entry) ? num_entry : max_num_entry;
91  }
92  max_num_entry *= 2; // 1 entry each for j, k coord
93  graph = Tpetra::createCrsGraph(map, max_num_entry);
94 
95  for (size_t i=0; i<pce_sz; ++i) {
96  const GlobalOrdinal row = i;
97  const size_t num_entry = cijk.num_entry(i);
98  const size_t entry_beg = cijk.entry_begin(i);
99  const size_t entry_end = entry_beg + num_entry;
100  for (size_t entry = entry_beg; entry < entry_end; ++entry) {
101  const GlobalOrdinal j = cijk.coord(entry,0);
102  const GlobalOrdinal k = cijk.coord(entry,1);
103  graph->insertGlobalIndices(row, arrayView(&j, 1));
104  graph->insertGlobalIndices(row, arrayView(&k, 1));
105  }
106  }
107  }
108  graph->fillComplete();
109  return graph;
110  }
111 
112  // Create a flattened graph for a graph from a matrix with the
113  // UQ::PCE scalar type
114  // If flat_domain_map and/or flat_range_map are null, they will be computed,
115  // otherwise they will be used as-is.
116  template <typename LocalOrdinal, typename GlobalOrdinal, typename Device,
117  typename CijkType>
118  Teuchos::RCP< Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,
119  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
120  create_flat_pce_graph(
121  const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& graph,
122  const CijkType& cijk,
123  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_domain_map,
124  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_range_map,
125  Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& cijk_graph,
126  const size_t matrix_pce_size) {
127  using Teuchos::ArrayView;
128  using Teuchos::ArrayRCP;
129  using Teuchos::Array;
130  using Teuchos::RCP;
131  using Teuchos::rcp;
132 
133  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
134  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
135  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
136 
137  const LocalOrdinal block_size = cijk.dimension();
138 
139  // Build domain map if necessary
140  if (flat_domain_map == Teuchos::null)
141  flat_domain_map = create_flat_map(*(graph.getDomainMap()), block_size);
142 
143  // Build range map if necessary
144  if (flat_range_map == Teuchos::null)
145  flat_range_map = create_flat_map(*(graph.getRangeMap()), block_size);
146 
147  // Build column map
148  RCP<const Map> flat_col_map =
149  create_flat_map(*(graph.getColMap()), block_size);
150 
151  // Build row map if necessary
152  // Check if range_map == row_map, then we can use flat_range_map
153  // as the flattened row map
154  RCP<const Map> flat_row_map;
155  if (graph.getRangeMap() == graph.getRowMap())
156  flat_row_map = flat_range_map;
157  else
158  flat_row_map = create_flat_map(*(graph.getRowMap()), block_size);
159 
160  // Build Cijk graph if necessary
161  if (cijk_graph == Teuchos::null)
162  cijk_graph = create_cijk_crs_graph<LocalOrdinal,GlobalOrdinal,Device>(
163  cijk,
164  flat_domain_map->getComm(),
165  matrix_pce_size);
166 
167  // Build flattened graph that is the Kronecker product of the given
168  // graph and cijk_graph
169 
170  // Loop over outer rows
171  ArrayView<const LocalOrdinal> outer_cols;
172  ArrayView<const LocalOrdinal> inner_cols;
173  size_t max_num_row_entries = graph.getNodeMaxNumRowEntries()*block_size;
174  Array<LocalOrdinal> flat_col_indices;
175  flat_col_indices.reserve(max_num_row_entries);
176  RCP<Graph> flat_graph = rcp(new Graph(flat_row_map, flat_col_map, max_num_row_entries));
177  const LocalOrdinal num_outer_rows = graph.getNodeNumRows();
178  for (LocalOrdinal outer_row=0; outer_row < num_outer_rows; outer_row++) {
179 
180  // Get outer columns for this outer row
181  graph.getLocalRowView(outer_row, outer_cols);
182  const LocalOrdinal num_outer_cols = outer_cols.size();
183 
184  // Loop over inner rows
185  for (LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
186 
187  // Compute flat row index
188  const LocalOrdinal flat_row = outer_row*block_size + inner_row;
189 
190  // Get inner columns for this inner row
191  cijk_graph->getLocalRowView(inner_row, inner_cols);
192  const LocalOrdinal num_inner_cols = inner_cols.size();
193 
194  // Create space to store all column indices for this flat row
195  flat_col_indices.resize(0);
196 
197  // Loop over outer cols
198  for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
199  ++outer_entry) {
200  const LocalOrdinal outer_col = outer_cols[outer_entry];
201 
202  // Loop over inner cols
203  for (LocalOrdinal inner_entry=0; inner_entry<num_inner_cols;
204  ++inner_entry) {
205  const LocalOrdinal inner_col = inner_cols[inner_entry];
206 
207  // Compute and store flat column index
208  const LocalOrdinal flat_col = outer_col*block_size + inner_col;
209  flat_col_indices.push_back(flat_col);
210  }
211 
212  }
213 
214  // Insert all indices for this flat row
215  flat_graph->insertLocalIndices(flat_row, flat_col_indices());
216 
217  }
218 
219  }
220  flat_graph->fillComplete(flat_domain_map, flat_range_map);
221 
222  return flat_graph;
223  }
224 
225  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
226  // returned vector is a view of the original
227  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
228  typename Device>
229  Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
230  LocalOrdinal,GlobalOrdinal,
231  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
232  create_flat_vector_view(
233  const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
234  LocalOrdinal,GlobalOrdinal,
235  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
236  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
237  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
238  using Teuchos::RCP;
239  using Teuchos::rcp;
240 
241  typedef typename Storage::value_type BaseScalar;
242  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
243  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
244  typedef typename FlatVector::dual_view_type flat_view_type;
245 
246  // Create flattenend view using special reshaping view assignment operator
247  flat_view_type flat_vals (vec.getLocalViewDevice(), vec.getLocalViewHost());
248  if (vec.need_sync_device ()) {
249  flat_vals.modify_host ();
250  }
251  else if (vec.need_sync_host ()) {
252  flat_vals.modify_device ();
253  }
254 
255  // Create flat vector
256  RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
257 
258  return flat_vec;
259  }
260 
261  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
262  // returned vector is a view of the original
263  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
264  typename Device>
265  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
266  LocalOrdinal,GlobalOrdinal,
267  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
268  create_flat_vector_view(
269  Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
270  LocalOrdinal,GlobalOrdinal,
271  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
272  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
273  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
274  using Teuchos::RCP;
275  using Teuchos::rcp;
276 
277  typedef typename Storage::value_type BaseScalar;
278  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
279  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
280  typedef typename FlatVector::dual_view_type flat_view_type;
281 
282  // Create flattenend view using special reshaping view assignment operator
283  flat_view_type flat_vals (vec.getLocalViewDevice(), vec.getLocalViewHost());
284  if (vec.need_sync_device ()) {
285  flat_vals.modify_host ();
286  }
287  else if (vec.need_sync_host ()) {
288  flat_vals.modify_device ();
289  }
290 
291  // Create flat vector
292  RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
293 
294  return flat_vec;
295  }
296 
297  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
298  // returned vector is a view of the original. This version creates the
299  // map if necessary
300  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
301  typename Device>
302  Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
303  LocalOrdinal,GlobalOrdinal,
304  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
305  create_flat_vector_view(
306  const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
307  LocalOrdinal,GlobalOrdinal,
308  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
309  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
310  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
311  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
312  if (flat_map == Teuchos::null) {
313  const LocalOrdinal pce_size =
314  Kokkos::dimension_scalar(vec.template getLocalView<Device>());
315  flat_map = create_flat_map(*(vec.getMap()), pce_size);
316  }
317  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
318  return create_flat_vector_view(vec, const_flat_map);
319  }
320 
321  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
322  // returned vector is a view of the original. This version creates the
323  // map if necessary
324  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
325  typename Device>
326  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
327  LocalOrdinal,GlobalOrdinal,
328  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
329  create_flat_vector_view(
330  Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
331  LocalOrdinal,GlobalOrdinal,
332  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
333  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
334  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
335  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
336  if (flat_map == Teuchos::null) {
337  const LocalOrdinal pce_size =
338  Kokkos::dimension_scalar(vec.template getLocalView<Device>());
339  flat_map = create_flat_map(*(vec.getMap()), pce_size);
340  }
341  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
342  return create_flat_vector_view(vec, const_flat_map);
343  }
344 
345  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
346  // returned vector is a view of the original
347  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
348  typename Device>
349  Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
350  LocalOrdinal,GlobalOrdinal,
351  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
352  create_flat_vector_view(
353  const Tpetra::Vector<Sacado::UQ::PCE<Storage>,
354  LocalOrdinal,GlobalOrdinal,
355  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec_const,
356  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
357  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
358  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
359  const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec_const;
360  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv = create_flat_vector_view(mv, flat_map);
361  return flat_mv->getVector(0);
362  }
363 
364  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
365  // returned vector is a view of the original. This version creates the
366  // map if necessary
367  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
368  typename Device>
369  Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
370  LocalOrdinal,GlobalOrdinal,
371  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
372  create_flat_vector_view(
373  const Tpetra::Vector<Sacado::UQ::PCE<Storage>,
374  LocalOrdinal,GlobalOrdinal,
375  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
376  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
377  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
378  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
379  if (flat_map == Teuchos::null) {
380  const LocalOrdinal pce_size =
381  Kokkos::dimension_scalar(vec.template getLocalView<Device>());
382  flat_map = create_flat_map(*(vec.getMap()), pce_size);
383  }
384  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
385  return create_flat_vector_view(vec, const_flat_map);
386  }
387 
388  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
389  // returned vector is a view of the original
390  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
391  typename Device>
392  Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
393  LocalOrdinal,GlobalOrdinal,
394  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
395  create_flat_vector_view(
396  Tpetra::Vector<Sacado::UQ::PCE<Storage>,
397  LocalOrdinal,GlobalOrdinal,
398  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
399  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
400  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
401  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
402  Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec;
403  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,LocalOrdinal,GlobalOrdinal,Node> > flat_mv = create_flat_vector_view(mv, flat_map);
404  return flat_mv->getVectorNonConst(0);
405  }
406 
407  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
408  // returned vector is a view of the original. This version creates the
409  // map if necessary
410  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
411  typename Device>
412  Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
413  LocalOrdinal,GlobalOrdinal,
414  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
415  create_flat_vector_view(
416  Tpetra::Vector<Sacado::UQ::PCE<Storage>,
417  LocalOrdinal,GlobalOrdinal,
418  Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
419  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,
420  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
421  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
422  if (flat_map == Teuchos::null) {
423  const LocalOrdinal pce_size =
424  Kokkos::dimension_scalar(vec.template getLocalView<Device>());
425  flat_map = create_flat_map(*(vec.getMap()), pce_size);
426  }
427  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
428  return create_flat_vector_view(vec, const_flat_map);
429  }
430 
431  // Create a flattened matrix by unrolling the UQ::PCE scalar type. The
432  // returned matrix is NOT a view of the original (and can't be)
433  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
434  typename Device, typename CijkType>
435  Teuchos::RCP< Tpetra::CrsMatrix<typename Storage::value_type,
436  LocalOrdinal,GlobalOrdinal,
437  Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
438  create_flat_matrix(
439  const Tpetra::CrsMatrix<Sacado::UQ::PCE<Storage>,
440  LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& mat,
441  const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_graph,
442  const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& cijk_graph,
443  const CijkType& cijk) {
444  using Teuchos::ArrayView;
445  using Teuchos::Array;
446  using Teuchos::RCP;
447  using Teuchos::rcp;
448 
449  typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
450  typedef Sacado::UQ::PCE<Storage> Scalar;
451  typedef typename Storage::value_type BaseScalar;
452  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
453 
454  const LocalOrdinal block_size = cijk.dimension();
455  const LocalOrdinal matrix_pce_size =
456  Kokkos::dimension_scalar(mat.getLocalMatrix().values);
457 
458  // Create flat matrix
459  RCP<FlatMatrix> flat_mat = rcp(new FlatMatrix(flat_graph));
460 
461  // Fill flat matrix
462  ArrayView<const Scalar> outer_values;
463  ArrayView<const LocalOrdinal> outer_cols;
464  ArrayView<const LocalOrdinal> inner_cols;
465  ArrayView<const LocalOrdinal> flat_cols;
466  Array<BaseScalar> flat_values;
467  flat_values.reserve(flat_graph->getNodeMaxNumRowEntries());
468  const LocalOrdinal num_outer_rows = mat.getNodeNumRows();
469  for (LocalOrdinal outer_row=0; outer_row < num_outer_rows; outer_row++) {
470 
471  // Get outer columns and values for this outer row
472  mat.getLocalRowView(outer_row, outer_cols, outer_values);
473  const LocalOrdinal num_outer_cols = outer_cols.size();
474 
475  // Loop over inner rows
476  for (LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
477 
478  // Compute flat row index
479  const LocalOrdinal flat_row = outer_row*block_size + inner_row;
480 
481  // Get cijk column indices for this row
482  cijk_graph->getLocalRowView(inner_row, inner_cols);
483  const LocalOrdinal num_inner_cols = inner_cols.size();
484 
485  // Create space to store all values for this flat row
486  const LocalOrdinal num_flat_indices = num_outer_cols*num_inner_cols;
487  //flat_values.resize(num_flat_indices);
488  flat_values.assign(num_flat_indices, BaseScalar(0));
489 
490  if (matrix_pce_size == 1) {
491  // Mean-based case
492 
493  // Loop over outer cols
494  for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
495  ++outer_entry) {
496 
497  // Extract mean PCE entry for each outer column
498  flat_values[outer_entry] = outer_values[outer_entry].coeff(0);
499 
500  }
501 
502  }
503  else {
504 
505  // Loop over cijk non-zeros for this inner row
506  const size_t num_entry = cijk.num_entry(inner_row);
507  const size_t entry_beg = cijk.entry_begin(inner_row);
508  const size_t entry_end = entry_beg + num_entry;
509  for (size_t entry = entry_beg; entry < entry_end; ++entry) {
510  const LocalOrdinal j = cijk.coord(entry,0);
511  const LocalOrdinal k = cijk.coord(entry,1);
512  const BaseScalar c = cijk.value(entry);
513 
514  // Find column offset for j
515  typedef typename ArrayView<const LocalOrdinal>::iterator iterator;
516  iterator ptr_j =
517  std::find(inner_cols.begin(), inner_cols.end(), j);
518  iterator ptr_k =
519  std::find(inner_cols.begin(), inner_cols.end(), k);
520  const LocalOrdinal j_offset = ptr_j - inner_cols.begin();
521  const LocalOrdinal k_offset = ptr_k - inner_cols.begin();
522 
523  // Loop over outer cols
524  for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
525  ++outer_entry) {
526 
527  // Add contributions for each outer column
528  flat_values[outer_entry*num_inner_cols + j_offset] +=
529  c*outer_values[outer_entry].coeff(k);
530  flat_values[outer_entry*num_inner_cols + k_offset] +=
531  c*outer_values[outer_entry].coeff(j);
532 
533  }
534 
535  }
536 
537  }
538 
539  // Set values in flat matrix
540  flat_graph->getLocalRowView(flat_row, flat_cols);
541  flat_mat->replaceLocalValues(flat_row, flat_cols, flat_values());
542 
543  }
544 
545  }
546  flat_mat->fillComplete(flat_graph->getDomainMap(),
547  flat_graph->getRangeMap());
548 
549  return flat_mat;
550  }
551 
552 #endif
553 
554 } // namespace Stokhos
555 
556 #endif // STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP
Storage
Stokhos::StandardStorage< int, double > Storage
Definition: linear2d_diffusion_scalar_types.hpp:47
Kokkos::dimension_scalar
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< View< T, P...> >::value, unsigned >::type dimension_scalar(const View< T, P...> &view)
Definition: KokkosExp_View_UQ_PCE_Contiguous.hpp:208
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_Tpetra_Utilities_MP_Vector.hpp
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
Kokkos::cijk
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Definition: KokkosExp_View_UQ_PCE_Contiguous.hpp:217
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
Stokhos_Sacado_Kokkos_UQ_PCE.hpp
Generated by   doxygen 1.8.5