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