10 #ifndef STOKHOS_TPETRA_UTILITIES_UQ_PCE_HPP
11 #define STOKHOS_TPETRA_UTILITIES_UQ_PCE_HPP
15 #include "Tpetra_Map.hpp"
16 #include "Tpetra_MultiVector.hpp"
17 #include "Tpetra_CrsGraph.hpp"
18 #include "Tpetra_CrsMatrix.hpp"
29 const size_t matrix_pce_size) {
31 using Teuchos::arrayView;
33 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
34 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
41 const size_t pce_sz =
cijk.dimension();
43 Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(pce_sz, comm);
45 if (matrix_pce_size == 1) {
46 graph = Tpetra::createCrsGraph(map, 1);
48 for (
size_t i=0; i<pce_sz; ++i) {
50 graph->insertGlobalIndices(row, arrayView(&row, 1));
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;
63 graph = Tpetra::createCrsGraph(map, max_num_entry);
65 for (
size_t i=0; i<pce_sz; ++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) {
73 graph->insertGlobalIndices(row, arrayView(&j, 1));
74 graph->insertGlobalIndices(row, arrayView(&k, 1));
78 graph->fillComplete();
90 const Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node >& graph,
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) {
102 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
103 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Graph;
109 flat_domain_map =
create_flat_map(*(graph.getDomainMap()), block_size);
123 if (graph.getRangeMap() == graph.getRowMap())
124 flat_row_map = flat_range_map;
130 cijk_graph = create_cijk_crs_graph<LocalOrdinal,GlobalOrdinal,Node>(
132 flat_domain_map->getComm(),
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;
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++) {
150 graph.getLocalRowView(outer_row, outer_cols);
154 for (
LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
157 const LocalOrdinal flat_row = outer_row*block_size + inner_row;
161 cijk_graph->getLocalRowView(inner_row, inner_cols);
165 flat_col_indices.
resize(0);
168 for (
LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
173 for (
LocalOrdinal inner_entry=0; inner_entry<num_inner_cols;
178 const LocalOrdinal flat_col = outer_col*block_size + inner_col;
185 flat_graph->insertLocalIndices(flat_row, flat_col_indices());
190 flat_graph->fillComplete(flat_domain_map, flat_range_map);
204 const Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
209 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
210 typedef typename FlatVector::dual_view_type::t_dev flat_view_type;
217 mv_type& vec_nc =
const_cast<mv_type&
>(vec);
220 flat_view_type flat_vals = vec_nc.getLocalViewDevice(Tpetra::Access::ReadWrite);
237 const Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
242 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
243 typedef typename FlatVector::dual_view_type::t_dev flat_view_type;
246 flat_view_type flat_vals = vec.getLocalViewDevice(Tpetra::Access::ReadWrite);
264 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
265 typedef typename Node::device_type Device;
285 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
286 typedef typename Node::device_type Device;
305 const Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
308 return flat_mv->getVector(0);
321 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
322 typedef typename Node::device_type Device;
341 const Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
344 return flat_mv->getVectorNonConst(0);
357 Teuchos::RCP<
const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node > >& flat_map) {
358 typedef typename Node::device_type Device;
371 typename Node,
typename CijkType>
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) {
387 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Matrix;
388 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
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;
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++) {
413 mat.getLocalRowView(outer_row, outer_cols, outer_values);
417 for (
LocalOrdinal inner_row=0; inner_row < block_size; inner_row++) {
420 const LocalOrdinal flat_row = outer_row*block_size + inner_row;
423 cijk_graph->getLocalRowView(inner_row, inner_cols);
426 Kokkos::Compat::getConstArrayView(inner_cols);
429 const LocalOrdinal num_flat_indices = num_outer_cols*num_inner_cols;
431 flat_values.
assign(num_flat_indices, BaseScalar(0));
433 if (matrix_pce_size == 1) {
437 for (
LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
441 flat_values[outer_entry] = outer_values[outer_entry].coeff(0);
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) {
455 const BaseScalar c =
cijk.value(entry);
460 std::find(inner_cols_av.
begin(), inner_cols_av.
end(),
j);
462 std::find(inner_cols_av.
begin(), inner_cols_av.
end(), k);
467 for (
LocalOrdinal outer_entry=0; outer_entry<num_outer_cols;
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);
483 flat_graph->getLocalRowView(flat_row, flat_cols);
484 flat_mat->replaceLocalValues(flat_row, Kokkos::Compat::getConstArrayView(flat_cols), flat_values());
489 flat_mat->fillComplete(flat_graph->getDomainMap(),
490 flat_graph->getRangeMap());
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
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)
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)