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 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef STOKHOS_TPETRA_UTILITIES_UQ_PCE_HPP
11 #define STOKHOS_TPETRA_UTILITIES_UQ_PCE_HPP
12 
15 #include "Tpetra_Map.hpp"
16 #include "Tpetra_MultiVector.hpp"
17 #include "Tpetra_CrsGraph.hpp"
18 #include "Tpetra_CrsMatrix.hpp"
19 
20 namespace Stokhos {
21 
22 
23  // Build a CRS graph from a sparse Cijk tensor
24  template <typename LocalOrdinal, typename GlobalOrdinal, typename Node,
25  typename CijkType>
27  create_cijk_crs_graph(const CijkType& cijk_dev,
28  const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
29  const size_t matrix_pce_size) {
30  using Teuchos::RCP;
31  using Teuchos::arrayView;
32 
33  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
34  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
35 
36  // Code below accesses cijk entries on the host, so make sure it is
37  // accessible there
38  auto cijk = create_mirror_view(cijk_dev);
39  deep_copy(cijk, cijk_dev);
40 
41  const size_t pce_sz = cijk.dimension();
43  Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(pce_sz, comm);
44  RCP<Graph> graph;
45  if (matrix_pce_size == 1) {
46  graph = Tpetra::createCrsGraph(map, 1);
47  // Mean-based case -- graph is diagonal
48  for (size_t i=0; i<pce_sz; ++i) {
49  const GlobalOrdinal row = i;
50  graph->insertGlobalIndices(row, arrayView(&row, 1));
51  }
52  }
53  else {
54  // General case
55 
56  // Get max num entries
57  size_t max_num_entry = 0;
58  for (size_t i=0; i<pce_sz; ++i) {
59  const size_t num_entry = cijk.num_entry(i);
60  max_num_entry = (num_entry > max_num_entry) ? num_entry : max_num_entry;
61  }
62  max_num_entry *= 2; // 1 entry each for j, k coord
63  graph = Tpetra::createCrsGraph(map, max_num_entry);
64 
65  for (size_t i=0; i<pce_sz; ++i) {
66  const GlobalOrdinal row = i;
67  const size_t num_entry = cijk.num_entry(i);
68  const size_t entry_beg = cijk.entry_begin(i);
69  const size_t entry_end = entry_beg + num_entry;
70  for (size_t entry = entry_beg; entry < entry_end; ++entry) {
71  const GlobalOrdinal j = cijk.coord(entry,0);
72  const GlobalOrdinal k = cijk.coord(entry,1);
73  graph->insertGlobalIndices(row, arrayView(&j, 1));
74  graph->insertGlobalIndices(row, arrayView(&k, 1));
75  }
76  }
77  }
78  graph->fillComplete();
79  return graph;
80  }
81 
82  // Create a flattened graph for a graph from a matrix with the
83  // UQ::PCE scalar type
84  // If flat_domain_map and/or flat_range_map are null, they will be computed,
85  // otherwise they will be used as-is.
86  template <typename LocalOrdinal, typename GlobalOrdinal, typename Node,
87  typename CijkType>
90  const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node >& graph,
91  const CijkType& cijk,
92  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_domain_map,
93  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_range_map,
94  Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node > >& cijk_graph,
95  const size_t matrix_pce_size) {
96  using Teuchos::ArrayView;
97  using Teuchos::ArrayRCP;
98  using Teuchos::Array;
99  using Teuchos::RCP;
100  using Teuchos::rcp;
101 
102  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
103  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
104 
105  const LocalOrdinal block_size = cijk.dimension();
106 
107  // Build domain map if necessary
108  if (flat_domain_map == Teuchos::null)
109  flat_domain_map = create_flat_map(*(graph.getDomainMap()), block_size);
110 
111  // Build range map if necessary
112  if (flat_range_map == Teuchos::null)
113  flat_range_map = create_flat_map(*(graph.getRangeMap()), block_size);
114 
115  // Build column map
116  RCP<const Map> flat_col_map =
117  create_flat_map(*(graph.getColMap()), block_size);
118 
119  // Build row map if necessary
120  // Check if range_map == row_map, then we can use flat_range_map
121  // as the flattened row map
122  RCP<const Map> flat_row_map;
123  if (graph.getRangeMap() == graph.getRowMap())
124  flat_row_map = flat_range_map;
125  else
126  flat_row_map = create_flat_map(*(graph.getRowMap()), block_size);
127 
128  // Build Cijk graph if necessary
129  if (cijk_graph == Teuchos::null)
130  cijk_graph = create_cijk_crs_graph<LocalOrdinal,GlobalOrdinal,Node>(
131  cijk,
132  flat_domain_map->getComm(),
133  matrix_pce_size);
134 
135  // Build flattened graph that is the Kronecker product of the given
136  // graph and cijk_graph
137 
138  // Loop over outer rows
139  typename Graph::local_inds_host_view_type outer_cols;
140  typename Graph::local_inds_host_view_type inner_cols;
141  size_t max_num_row_entries = graph.getLocalMaxNumRowEntries()*block_size;
142  Array<LocalOrdinal> flat_col_indices;
143  flat_col_indices.reserve(max_num_row_entries);
144  RCP<Graph> flat_graph = rcp(new Graph(flat_row_map, flat_col_map, max_num_row_entries));
145  const LocalOrdinal num_outer_rows = graph.getLocalNumRows();
146  for (LocalOrdinal outer_row=0; outer_row < num_outer_rows; outer_row++) {
147 
148  // Get outer columns for this outer row
149  Kokkos::fence();
150  graph.getLocalRowView(outer_row, outer_cols);
151  const LocalOrdinal num_outer_cols = outer_cols.size();
152 
153  // Loop over inner rows
154  for (LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
155 
156  // Compute flat row index
157  const LocalOrdinal flat_row = outer_row*block_size + inner_row;
158 
159  // Get inner columns for this inner row
160  Kokkos::fence();
161  cijk_graph->getLocalRowView(inner_row, inner_cols);
162  const LocalOrdinal num_inner_cols = inner_cols.size();
163 
164  // Create space to store all column indices for this flat row
165  flat_col_indices.resize(0);
166 
167  // Loop over outer cols
168  for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
169  ++outer_entry) {
170  const LocalOrdinal outer_col = outer_cols[outer_entry];
171 
172  // Loop over inner cols
173  for (LocalOrdinal inner_entry=0; inner_entry<num_inner_cols;
174  ++inner_entry) {
175  const LocalOrdinal inner_col = inner_cols[inner_entry];
176 
177  // Compute and store flat column index
178  const LocalOrdinal flat_col = outer_col*block_size + inner_col;
179  flat_col_indices.push_back(flat_col);
180  }
181 
182  }
183 
184  // Insert all indices for this flat row
185  flat_graph->insertLocalIndices(flat_row, flat_col_indices());
186 
187  }
188 
189  }
190  flat_graph->fillComplete(flat_domain_map, flat_range_map);
191 
192  return flat_graph;
193  }
194 
195  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
196  // returned vector is a view of the original
197  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
198  typename Node>
199  Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
202  const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
204  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
205  using Teuchos::RCP;
206  using Teuchos::rcp;
207 
208  typedef typename Storage::value_type BaseScalar;
209  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
210  typedef typename FlatVector::dual_view_type::t_dev flat_view_type;
211 
212  // Have to do a nasty const-cast because getLocalViewDevice(ReadWrite) is a
213  // non-const method, yet getLocalViewDevice(ReadOnly) returns a const-view
214  // (i.e., with a constant scalar type), and there is no way to make a
215  // MultiVector out of it!
216  typedef Tpetra::MultiVector<Sacado::UQ::PCE<Storage>, LocalOrdinal,GlobalOrdinal, Node > mv_type;
217  mv_type& vec_nc = const_cast<mv_type&>(vec);
218 
219  // Create flattenend view using special reshaping view assignment operator
220  flat_view_type flat_vals = vec_nc.getLocalViewDevice(Tpetra::Access::ReadWrite);
221 
222  // Create flat vector
223  RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
224 
225  return flat_vec;
226  }
227 
228  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
229  // returned vector is a view of the original
230  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
231  typename Node>
232  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
235  Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
237  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
238  using Teuchos::RCP;
239  using Teuchos::rcp;
240 
241  typedef typename Storage::value_type BaseScalar;
242  typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
243  typedef typename FlatVector::dual_view_type::t_dev flat_view_type;
244 
245  // Create flattenend view using special reshaping view assignment operator
246  flat_view_type flat_vals = vec.getLocalViewDevice(Tpetra::Access::ReadWrite);
247 
248  // Create flat vector
249  RCP<FlatVector> flat_vec = rcp(new FlatVector(flat_map, flat_vals));
250 
251  return flat_vec;
252  }
253 
254  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
255  // returned vector is a view of the original. This version creates the
256  // map if necessary
257  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
258  typename Node>
259  Teuchos::RCP< const Tpetra::MultiVector<typename Storage::value_type,
262  const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
264  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
265  typedef typename Node::device_type Device;
266  if (flat_map == Teuchos::null) {
267  const LocalOrdinal pce_size =
268  Kokkos::dimension_scalar(vec.template getLocalView<Device>(Tpetra::Access::ReadOnly));
269  flat_map = create_flat_map(*(vec.getMap()), pce_size);
270  }
272  return create_flat_vector_view(vec, const_flat_map);
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 Node>
280  Teuchos::RCP< Tpetra::MultiVector<typename Storage::value_type,
283  Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
285  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
286  typedef typename Node::device_type Device;
287  if (flat_map == Teuchos::null) {
288  const LocalOrdinal pce_size =
289  Kokkos::dimension_scalar(vec.template getLocalView<Device>(Tpetra::Access::ReadOnly));
290  flat_map = create_flat_map(*(vec.getMap()), pce_size);
291  }
293  return create_flat_vector_view(vec, const_flat_map);
294  }
295 
296  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
297  // returned vector is a view of the original
298  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
299  typename Node>
300  Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
303  const Tpetra::Vector<Sacado::UQ::PCE<Storage>,
304  LocalOrdinal,GlobalOrdinal,Node >& vec_const,
305  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
306  const Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec_const;
308  return flat_mv->getVector(0);
309  }
310 
311  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
312  // returned vector is a view of the original. This version creates the
313  // map if necessary
314  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
315  typename Node>
316  Teuchos::RCP< const Tpetra::Vector<typename Storage::value_type,
319  const Tpetra::Vector<Sacado::UQ::PCE<Storage>,
321  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
322  typedef typename Node::device_type Device;
323  if (flat_map == Teuchos::null) {
324  const LocalOrdinal pce_size =
325  Kokkos::dimension_scalar(vec.template getLocalView<Device>(Tpetra::Access::ReadOnly));
326  flat_map = create_flat_map(*(vec.getMap()), pce_size);
327  }
329  return create_flat_vector_view(vec, const_flat_map);
330  }
331 
332  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
333  // returned vector is a view of the original
334  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
335  typename Node>
336  Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
339  Tpetra::Vector<Sacado::UQ::PCE<Storage>,
341  const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
342  Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,LocalOrdinal,GlobalOrdinal,Node>& mv = vec;
344  return flat_mv->getVectorNonConst(0);
345  }
346 
347  // Create a flattened vector by unrolling the UQ::PCE scalar type. The
348  // returned vector is a view of the original. This version creates the
349  // map if necessary
350  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
351  typename Node>
352  Teuchos::RCP< Tpetra::Vector<typename Storage::value_type,
355  Tpetra::Vector<Sacado::UQ::PCE<Storage>,
357  Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
358  typedef typename Node::device_type Device;
359  if (flat_map == Teuchos::null) {
360  const LocalOrdinal pce_size =
361  Kokkos::dimension_scalar(vec.template getLocalView<Device>(Tpetra::Access::ReadOnly));
362  flat_map = create_flat_map(*(vec.getMap()), pce_size);
363  }
365  return create_flat_vector_view(vec, const_flat_map);
366  }
367 
368  // Create a flattened matrix by unrolling the UQ::PCE scalar type. The
369  // returned matrix is NOT a view of the original (and can't be)
370  template <typename Storage, typename LocalOrdinal, typename GlobalOrdinal,
371  typename Node, typename CijkType>
372  Teuchos::RCP< Tpetra::CrsMatrix<typename Storage::value_type,
375  const Tpetra::CrsMatrix<Sacado::UQ::PCE<Storage>,
377  const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node > >& flat_graph,
378  const Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node > >& cijk_graph,
379  const CijkType& cijk_dev) {
380  using Teuchos::ArrayView;
381  using Teuchos::Array;
382  using Teuchos::RCP;
383  using Teuchos::rcp;
384 
386  typedef typename Storage::value_type BaseScalar;
387  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Matrix;
388  typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
389 
390  // Code below accesses cijk entries on the host, so make sure it is
391  // accessible there
392  auto cijk = create_mirror_view(cijk_dev);
393  deep_copy(cijk, cijk_dev);
394 
395  const LocalOrdinal block_size = cijk.dimension();
396  const LocalOrdinal matrix_pce_size =
397  Kokkos::dimension_scalar(mat.getLocalMatrixDevice().values);
398 
399  // Create flat matrix
400  RCP<FlatMatrix> flat_mat = rcp(new FlatMatrix(flat_graph));
401 
402  // Fill flat matrix
403  typename Matrix::values_host_view_type outer_values;
404  typename Matrix::local_inds_host_view_type outer_cols;
405  typename Matrix::local_inds_host_view_type inner_cols;
406  typename Matrix::local_inds_host_view_type flat_cols;
407  Array<BaseScalar> flat_values;
408  flat_values.reserve(flat_graph->getLocalMaxNumRowEntries());
409  const LocalOrdinal num_outer_rows = mat.getLocalNumRows();
410  for (LocalOrdinal outer_row=0; outer_row < num_outer_rows; outer_row++) {
411 
412  // Get outer columns and values for this outer row
413  mat.getLocalRowView(outer_row, outer_cols, outer_values);
414  const LocalOrdinal num_outer_cols = outer_cols.size();
415 
416  // Loop over inner rows
417  for (LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
418 
419  // Compute flat row index
420  const LocalOrdinal flat_row = outer_row*block_size + inner_row;
421 
422  // Get cijk column indices for this row
423  cijk_graph->getLocalRowView(inner_row, inner_cols);
424  const LocalOrdinal num_inner_cols = inner_cols.size();
425  ArrayView<const LocalOrdinal> inner_cols_av =
426  Kokkos::Compat::getConstArrayView(inner_cols);
427 
428  // Create space to store all values for this flat row
429  const LocalOrdinal num_flat_indices = num_outer_cols*num_inner_cols;
430  //flat_values.resize(num_flat_indices);
431  flat_values.assign(num_flat_indices, BaseScalar(0));
432 
433  if (matrix_pce_size == 1) {
434  // Mean-based case
435 
436  // Loop over outer cols
437  for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
438  ++outer_entry) {
439 
440  // Extract mean PCE entry for each outer column
441  flat_values[outer_entry] = outer_values[outer_entry].coeff(0);
442 
443  }
444 
445  }
446  else {
447 
448  // Loop over cijk non-zeros for this inner row
449  const size_t num_entry = cijk.num_entry(inner_row);
450  const size_t entry_beg = cijk.entry_begin(inner_row);
451  const size_t entry_end = entry_beg + num_entry;
452  for (size_t entry = entry_beg; entry < entry_end; ++entry) {
453  const LocalOrdinal j = cijk.coord(entry,0);
454  const LocalOrdinal k = cijk.coord(entry,1);
455  const BaseScalar c = cijk.value(entry);
456 
457  // Find column offset for j
458  typedef typename ArrayView<const LocalOrdinal>::iterator iterator;
459  iterator ptr_j =
460  std::find(inner_cols_av.begin(), inner_cols_av.end(), j);
461  iterator ptr_k =
462  std::find(inner_cols_av.begin(), inner_cols_av.end(), k);
463  const LocalOrdinal j_offset = ptr_j - inner_cols_av.begin();
464  const LocalOrdinal k_offset = ptr_k - inner_cols_av.begin();
465 
466  // Loop over outer cols
467  for (LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
468  ++outer_entry) {
469 
470  // Add contributions for each outer column
471  flat_values[outer_entry*num_inner_cols + j_offset] +=
472  c*outer_values[outer_entry].coeff(k);
473  flat_values[outer_entry*num_inner_cols + k_offset] +=
474  c*outer_values[outer_entry].coeff(j);
475 
476  }
477 
478  }
479 
480  }
481 
482  // Set values in flat matrix
483  flat_graph->getLocalRowView(flat_row, flat_cols);
484  flat_mat->replaceLocalValues(flat_row, Kokkos::Compat::getConstArrayView(flat_cols), flat_values());
485 
486  }
487 
488  }
489  flat_mat->fillComplete(flat_graph->getDomainMap(),
490  flat_graph->getRangeMap());
491 
492  return flat_mat;
493  }
494 
495 
496 } // namespace Stokhos
497 
498 #endif // STOKHOS_TPETRA_UTILITIES_MP_VECTOR_HPP
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > create_cijk_crs_graph(const CijkType &cijk_dev, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const size_t matrix_pce_size)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Node > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &cijk_graph, const CijkType &cijk_dev)
void reserve(size_type n)
Stokhos::StandardStorage< int, double > Storage
iterator begin() const
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
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_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
CrsProductTensor< ValueType, Device, Memory >::HostMirror create_mirror_view(const CrsProductTensor< ValueType, Device, Memory > &src)
void resize(size_type new_size, const value_type &x=value_type())
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)
void push_back(const value_type &x)
void deep_copy(const CrsProductTensor< ValueType, DstDevice, DstMemory > &dst, const CrsProductTensor< ValueType, SrcDevice, SrcMemory > &src)
iterator end() const
Tpetra::KokkosClassic::DefaultNode::DefaultNodeType Node
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > create_flat_pce_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, const CijkType &cijk, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_range_map, Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &cijk_graph, const size_t matrix_pce_size)
Teuchos::RCP< Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > create_flat_map(const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map, const LocalOrdinal block_size)
void assign(size_type n, const value_type &val)