42 #ifndef STOKHOS_TPETRA_UTILITIES_UQ_PCE_HPP 43 #define STOKHOS_TPETRA_UTILITIES_UQ_PCE_HPP 47 #include "Tpetra_Map.hpp" 48 #include "Tpetra_MultiVector.hpp" 49 #include "Tpetra_CrsGraph.hpp" 50 #include "Tpetra_CrsMatrix.hpp" 54 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR) 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) {
66 using Teuchos::arrayView;
68 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
69 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
70 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
72 const size_t pce_sz =
cijk.dimension();
74 Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal>(pce_sz, comm, node);
75 RCP<Graph> graph = Tpetra::createCrsGraph(map);
76 if (matrix_pce_size == 1) {
78 for (
size_t i=0; i<pce_sz; ++i) {
80 graph->insertGlobalIndices(row, arrayView(&row, 1));
85 for (
size_t i=0; i<pce_sz; ++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) {
93 graph->insertGlobalIndices(row, arrayView(&
j, 1));
94 graph->insertGlobalIndices(row, arrayView(&k, 1));
98 graph->fillComplete();
109 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
110 create_flat_pce_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;
123 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
124 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
125 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
130 if (flat_domain_map == Teuchos::null)
131 flat_domain_map =
create_flat_map(*(graph.getDomainMap()), block_size);
134 if (flat_range_map == Teuchos::null)
138 RCP<const Map> flat_col_map =
144 RCP<const Map> flat_row_map;
145 if (graph.getRangeMap() == graph.getRowMap())
146 flat_row_map = flat_range_map;
151 if (cijk_graph == Teuchos::null)
152 cijk_graph = create_cijk_crs_graph<LocalOrdinal,GlobalOrdinal>(
154 flat_domain_map->getComm(),
155 flat_domain_map->getNode(),
160 RCP<Graph> flat_graph = rcp(
new Graph(flat_row_map, flat_col_map, 0));
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++) {
171 graph.getLocalRowView(outer_row, outer_cols);
175 for (
LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
178 const LocalOrdinal flat_row = outer_row*block_size + inner_row;
181 cijk_graph->getLocalRowView(inner_row, inner_cols);
185 flat_col_indices.resize(0);
188 for (
LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
193 for (
LocalOrdinal inner_entry=0; inner_entry<num_inner_cols;
198 const LocalOrdinal flat_col = outer_col*block_size + inner_col;
199 flat_col_indices.push_back(flat_col);
205 flat_graph->insertLocalIndices(flat_row, flat_col_indices());
210 flat_graph->fillComplete(flat_domain_map, flat_range_map);
221 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
225 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
227 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
232 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
234 typedef typename FlatVector::dual_view_type flat_view_type;
237 flat_view_type flat_vals = vec.getDualView();
240 RCP<FlatVector> flat_vec = rcp(
new FlatVector(flat_map, flat_vals));
251 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
255 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
257 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
262 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
264 typedef typename FlatVector::dual_view_type flat_view_type;
267 flat_view_type flat_vals = vec.getDualView();
270 RCP<FlatVector> flat_vec = rcp(
new FlatVector(flat_map, flat_vals));
282 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
286 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
288 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
289 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
290 if (flat_map == Teuchos::null) {
295 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
306 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
310 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
312 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
313 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
314 if (flat_map == Teuchos::null) {
319 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
329 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
333 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec_const,
335 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
336 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
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);
349 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
353 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
355 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
356 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
357 if (flat_map == Teuchos::null) {
362 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
372 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
376 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
378 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
379 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
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);
392 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
396 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >& vec,
398 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >& flat_map) {
399 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
400 if (flat_map == Teuchos::null) {
405 const Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > const_flat_map = flat_map;
412 typename Device,
typename CijkType>
415 Kokkos::Compat::KokkosDeviceWrapperNode<Device> > >
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;
427 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device>
Node;
430 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
437 RCP<FlatMatrix> flat_mat = rcp(
new FlatMatrix(flat_graph));
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++) {
450 mat.getLocalRowView(outer_row, outer_cols, outer_values);
454 for (
LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
457 const LocalOrdinal flat_row = outer_row*block_size + inner_row;
460 cijk_graph->getLocalRowView(inner_row, inner_cols);
464 const LocalOrdinal num_flat_indices = num_outer_cols*num_inner_cols;
466 flat_values.assign(num_flat_indices, BaseScalar(0));
468 if (matrix_pce_size == 1) {
472 for (
LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
476 flat_values[outer_entry] = outer_values[outer_entry].coeff(0);
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) {
490 const BaseScalar c =
cijk.value(entry);
493 typedef typename ArrayView<const LocalOrdinal>::iterator iterator;
495 std::find(inner_cols.begin(), inner_cols.end(),
j);
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();
502 for (
LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
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);
518 flat_graph->getLocalRowView(flat_row, flat_cols);
519 flat_mat->replaceLocalValues(flat_row, flat_cols, flat_values());
524 flat_mat->fillComplete(flat_graph->getDomainMap(),
525 flat_graph->getRangeMap());
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)
Teuchos::RCP< Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > create_flat_map(const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map, const LocalOrdinal block_size)